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First  principles  fully  nonlocal  pseudopotentials  are  constructed 
for  the  alkali  metals  Li,  Na,  and  K.  The  orthogonal ization  hole 
contribution  to  the  pseudopotential  is  treated  exactly  and  comparison 
with  approximate  treatments  show  significant  differences.  The  pseudo- 
potentials are  used  to  calculate  various  solid  and  liquid  properties 
of  these  alkalis.  The  phonon  spectra  and  elastic  shear  constants  are 
calculated  for  the  solid  metal  and  are  in  good  agreement  with  experi- 
ment. They  are  generally  within  20%  of  experiment  and  in  some  cases 
much  better  agreement  is  obtained.  The  liquid  structure  factor  curve, 

S (q) , is  calculated  using  a Monte  Carlo  technique  and  the  agreement 
with  experiment  is  excellent.  Electrical  resistivities  of  the  liquid 
alkali  metals  are  also  calculated  using  the  pseudopotential  and 
computed  liquid  structure.  The  electronic  transport  properties  are 
calculated  for  the  first  time,  using  the  same  ion  potential  to  determine 
the  theoretical  liquid  structure  and  to  describe  the  electron  scattering 
potential  in  the  solid.  The  results  are  in  reasonable  agreement  with 
experiment  considering  the  first  principles  nature  of  the  calculation. 

Finally,  a pseudo-atom  model  is  suggested  as  a more  appropriate 
means  of  describing  the  core-electron  states  used  in  constructing  the 
pseudopotential . 
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I.  INTRODUCTION  AND  STATEMENT  OF  OBJECTIVES 

In  this  thesis,  pseudopotential  calculations  are  done  for  the 
simple  (i.e.,  no  d-electrons)  alkali  metals  Li,  Na,  and  K.  These  metals 
were  chosen  because  they  are  free-electron-like  in  their  behavior  and 
have  a small  well  defined  core.  They  are,  therefore,  well  suited  for 
investigation  by  the  pseudopotential  method.  Solid  and  liquid  metal 
properties  are  calculated  for  each  of  the  aforementioned  alkali  metals. 
The  solid  properties  calculated  are  the  phonon  spectra,  the  pair 
potentials,  and  the  elastic  shear  constants.  The  liquid  metal  proper- 
ties calculated  include  the  effective  ion-ion  potential,  the  radial 
distribution  function,  the  static  structure  factor,  and  the  electrical 
resistivity  of  the  metal.  Both  the  solid  and  liquid  metal  results  are 
in  good  agreement  with  experimental  evidence. 

The  pseudopotentials  used  in  these  calculations  are  distinctive 
in  a number  of  ways.  First,  all  calculations  are  fully  nonlocal  through- 
out all  stages  of  the  computation.  The  only  exception  is  the  calculation 
of  the  electrical  resistivity  for  the  liquid  metal  which  uses  the  Ziman 
(1961)  local  theory  discussed  in  Appendix  E.  The  nonlocality  of  the 
pseudopotential  means  that  its  matrix  element  <k+q|w|k>  depends  on  the 
magnitude  of  k as  well  as  the  angle  between  k and  k+q,  and  is  not 
restricted  to  the  so-called  on  Fermi  sphere  approximation  (i.e., 

|k|  = k^  ).  Use  of  the  nonlocal  form  of  the  pseudopotential  makes  the 
computer  calculations  considerably  more  difficult,  but  computed  results 
can  be  quite  different  if  the  nonlocality  of  the  pseudopotential  is 


ignored.  Secondly,  it  should  be  stressed  that  our  pseudopotentials 


are  first  principles  which  ate  distinctly  different  than  model  ot 
pat arnate r i zed  potentials.  The  following  criteria  are  implied  by  the 
use  of  the  words  ’first  principles.’  A pseudo|<ot  ent  1 a 1 is  t i rst  prin- 
ciples it  the  localized  core-electron  states  arc  explicitly  represented 
by  wave  functions  that  are  solutions  of  a Schrodinger  equation.  These 
core  states  are  used  t.o  construct  the  ionic  charge  density  and, 
subsequently,  this  charge  density  determines  the  Coulombic  potential 
of  the  ion  in  the  crystal.  In  addition,  the  crystalline  potential  is 
represented  by  an  expression  that  does  not  make  use  of  arbitrarily 
adjustable  parameters.  The  only  experimental  information  used  in  our 
calculations,  other  than  fundamental  physical  constants,  is  the  lattice 
structure,  lattice  constant,  and  the  atomic  number  of  the  atoms  composing 
the  metal. 

The  basic  objective  of  this  thesis  is  t o study  the  electron-electron 
and  electron-ion  interactions  which  are  responsible  for  determining  the 
physical  properties  of  the  alkali  metals  Li,  Na,  and  K.  This  is  done 
by  determining  the  sensitivity  of  the  lattice  dynamical  properties  t o 
each  of  the  various  contributions  to  the  crystalline  potential.  As  a 
result  ot  these  studies  certain  contributions  to  the  pseudopot ent i a l 
ate  treated  more  rigorously  than  in  previous  calculat ions . Speci t ical ly , 
flu'  charge  density  and  Coulomb  potential  of  the  ort hogona 1 i zat i on  hole 
(i'll),  as  well  as  its  interaction  with  the  core-electron  density  are 
treated  exact  ly  for  the  first  time  for  Li,  Na,  and  K.  The  results 
using  an  ’exact’  i'll  density  show  that  approximate  treatments  ot  this 
particular  contribution  to  the  conduct  ion- elect  ton  potential  are 
inadequate  and  that  an  exact  treatment  of  the  Oil  produces  significant 


changes  in  tlu*  calculated  phonon  spectra.  Physically,  the  OH  density, 
n*-  **  U ) , is  a part  of  the  conduction-electron  charge  density  which  is 
obtained  tiom  the  absolute  square  of  the  conduct  ion-electron  wave 
function.  Also,  a thorough  analysis  is  made  of  the  conduction-core 
exchange  interaction  and  it  is  found  that  the  Lindgren  11971)  exchange 
potential  yields  better  results  in  most  cases,  than  the  Kohn-Sham 
(19bS)  or  Slatet  (19‘.'1)  potentials.  Finally,  the  electronic  transport 
properties  of  the  liquid  alk.nl  metals  are  calculated  for  the  first  time 
using  the  same  ion  potential  to  determine  the  theoretical  liquid 
structure  and  to  describe  the  electron  scattering  potential  in  the 


liquid. 
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11.  HISTORICAL  NOTH 

Of  central  importance  to  pseudopotential  theory  are  the  Ortho- 
gonalized  Plane  Waves  (Ol’W)  first  introduced  by  Herring  (1940).  These 
Ol'W's  are  used  to  represent  the  conduction-electron  states  in  a crystal. 
The  OFW  representation  for  the  conduction  eigenstates  has  the  advantage 
that  it  converges  much  faster  than  a simple  plane  wave  representation. 
Also,  it  was  noted  by  Phillips  and  Kleinman  (1959)  that  the  orthogona 1 iza- 
tion  of  the  plane  waves  to  the  core-electron  states  effectively  adds  a 
repulsive  term  to  the  crystalline  potential.  The  attractive  crystalline 
potential  plus  the  repulsive  orthogonal izat ion  term  have  been  collectively 
labelled  as  the  'pseudopotential . 1 This  potential  is  weaker  than  the 
crystalline  potential  which  suggests  that  one  can  use  perturbation  theory 
to  solve  for  the  conduction-electron  eigenstates.  Harrison  (1963,1963, 
1964)  was  the  first  to  apply  the  OPW  pseudopotential  theory  to  the 
calculation  of  metallic  properties.  His  calculations  are  noteworthy 
because  he  developed  his  psoudopotentials  from  first  principles.  This 
means  that  the  potential  is  developed  from  the  wave  functions  and 
eigenvalues  used  to  represent  the  localized  core-electron  states.  Any 
pseudopotent ial  in  which  the  crystalline  jiotential  depends  explicitly 
upon  the  core-electron  states  will  hereafter  be  referred  to  as  a 
Harrison's  First  Principles  (HFP)  pseudopotent ial . 

It  is  instructive  to  compare  the  HFP  pseudopotential  with  t tie 
'model*  pseudopotential  method,  which  does  not  construct  the  potential 


of  an  ion  in  the  crystal  from  the  explicit  core-electron  eigenstates. 


5 


Rather,  as  the  name  implies,  a model  potential  is  constructed  to 
represent  the  ionic  cores  in  the  metal . In  both  the  HFP  and  model 
pseudopotentials,  the  metal  is  visualized  as  an  array  of  ions  immersed 
in  a compensating  conduction-electron  sea.  Heine  and  Abrenkov  (1965), 
Animalu  (1965) , and  Ashcroft  (1966)  were  responsible  for  much  of  the 
early  development  of  model  pseudopotentials.  Heine  and  Weire  (1970) 
have  written  an  excellent  review  article  on  the  model  pseudopotential 
method,  so  we  will  only  briefly  touch  on  the  relevant  points  of  this 
approach.  In  such  calculations, the  potential  of  an  ion  is  represented 
by  a model  (i.e.,  an  analytic  expression  such  as  a square  well  potential) 
with  several  adjustable  parameters  that  are  used  to  fit  the  model  to 
known  experimental  facts.  It  is  hoped  that  once  the  parameters  have 
been  adjusted  to  one  measurement,  then  the  pseudopotential  can  be  used 
to  predict  other  atomic  and  electronic  properties  of  the  metal  under 
study. 

There  are  certain  model  potentials  that  are  quite  sophisticated 
in  the  way  that  the  parameters  are  chosen.  In  particular,  Dagens, 
Rasolt,  and  Taylor  (1975)  have  developed  a clever  approach  based  on 
the  charge  density  induced  by  an  isolated  ion  in  an  electron  gas.  They 
adjusted  the  model  potential  parameters  so  that  the  induced  charge 
density  calculated  from  first  order  perturbation  theory  was  identical 
with  an  exact  nonlinear  calculation  (Dagens,  1972).  This  is  certainly 
a more  fundamental  approach  than  adjusting  the  model  potential  to  the 
phonon  spectra,  elastic  constants,  or  other  gross  experimental  features. 
However,  irrespective  of  the  ingenuity  used  to  adjust  the  model  parame- 
ters, one  is  still  working  with  a model  potential  and  consequently  can 
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not  directly  obtain  fundamental  information  about  the  interactions  in 
the  crystal.  The  HFP  potential  by  contrast  constructs  the  crystal 
potential  from  wave  functions  that  are  solutions  of  the  Schrodinger 
equation,  and  thus  provides  a better  means  of  studying  the  crystal 
potential  and  the  interactions  that  contribute  to  the  crystal  potential 
than  does  the  model  pseudopotential  method. 

Though  the  HFP  pseudopotential  method  can,  in  principle,  be  used 
to  make  systematic  studies  of  the  interactions  existing  within  a 
crystal,  there  have  been  relatively  few  such  studies  done  because  of 
the  cost  and  complexity  of  the  calculations.  The  difficulty  of  the 
HFP  calculations  has  limited  the  practical  applications  of  this  tech- 
nique so  that  the  great  majority  of  pseudopotential  calculations  use 
the  model  potential  approach.  Since  the  basic  objective  of  this 
research  is  to  try  to  understand  some  fundamental  interactions  in  the 
solid  and  liquid  states  of  metals,  the  HFP  approach  has  been  chosen  to 
systematically  study  the  atomic  and  electronic  properties  of  the  alkali 
metals  Li,  Na,  and  K.  The  heavier  alkalis  are  not  included  in  this 
study  because  the  presence  of  d-electrons  in  Rb  and  Cs  may  require  the 
use  of  the  more  complicated  Generalized  HFP  (GHFP)  theory  which  includes 
the  s-d  'hybridization'  in  the  pseudopotential  (Harrison,  1969;  Moriarity, 
1970) . The  alkalis  are  ideally  suited  for  study  by  the  pseudopotential 
method  because  of  their  nearly  free  electron-like  behavior  and  the  fact 
that  there  is  virtually  no  nearest-neighbor  core  overlap  in  the  crystal. 
King  and  Cutler  (1973)  have  done  studies  for  Li  that  demonstrated  the 
importance  of  the  conduction-core  electron  exchange  interaction  in  phonon 
spectra  calculations.  More  recently,  Hafner  (1974) , and  Hafner  and 
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I 

Schmuck  (1974)  have  done  studies  which  also  conclude  that  the  conduction- 

core  exchange  interaction  is  an  important  component  of  the  crystal 

potential.  Hafner's  work,  however,  tended  to  use  parameterized  forms 

of  the  exchange  potential.  By  contrast,  the  present  work  represents  a 

systematic  study  of  the  exchange  potential  in  Li,  Na,  and  K using 

unparameterized  versions  of  approximations  to  the  conduction-core 

exchange  interaction.  In  addition,  for  the  first  time,  certain 

conduction-electron  contributions  to  the  crystal  potential  are  treated 

without  approximations;  this  exact  treatment  of  the  conduction- 

electron  potential  has  not  been  done  before  for  the  alkalis.  Spe- 

OH  "*• 

cifically,  the  Orthogonalization  Hole  density  n (r)  is  treated 
exactly  and  is  found  to  produce  significantly  different  results  from 
previous  approximate  treatments  of  n^H(r). 
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III.  THE  LINEAR  OPTIMIZED  PSEUDOPOTENTIAL 


Pseudopotential  Matrix  Elements 


In  this  section,  a derivation  is  given  of  the  pseudopotential  matrix 

elements  that  are  used  in  HFP  calculations  of  the  conduction-electron 

eigenvalues  and  wave  functions.  The  pseudopotential,  denoted  by  W(r), 

c 

will  be  expressed  in  terms  of  the  crystal  potential  V (r) , the  core- 
electron  eigenstates  (/^(r)  > and  their  corresponding  eigenvalues  E 
where  n and  l are  the  principle  and  angular  quantum  numbers  respectively. 
The  specific  forms  of  V (r) , ^n£<r)  and  En£  will  be  discussed  in  later 
chapters . 


In  HFP  pseudopotential  theory,  the  conduction-electron  states  are 


described  by  linear  combinations  of  OPW's,  where  a single  OPW  of  wave 
->  -> 

number  k + q can  be  written  as  (assuming  no  core  electron  overlaps) 

| OPW>  = (1  - P)  |k  + q>  , (1) 

and  P = l |n£><n£|  is  a projection  operator,  which  when  operating  on  the 

ni  i-  -+ 

plane  wave  state  |k  + q>,  extracts  the  components  of  the  plane  wave 
that  are  not  orthogonal  to  the  core-electron  states  |n£>.  Hence,  the 
OPW  in  Equation  (1)  represents  a plane  wave  that  has  been  made  orthogonal 
to  the  core-electron  states.  The  general  expression  for  a conduction- 
electron  state  of  wave  vector  k is  a linear  sum  of  the  OPW's  shown  in 


Equation  (1) 


i|£(r)  = l a (k)  (1  - P)  |k  + q> 


The  expansion  coefficients  a (k)  in  Equation  (2)  are  obtained  from 

q 

perturbation  theory. 


J 
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The  coefficients  are  determined  by  the  requirement  that  (r)  is 

k 

an  eigenstate  of  the  crystal  Hamiltonian  H = T + V . 


H^(r)  = Ej£(r)  (3) 

k k k 

Combining  Equations  (2)  and  (3),  and  with  a certain  amount  of  algebra 
(Harrison,  1966)  one  can  show  that  Equation  (3)  may  be  reexpressed  as 


(T  + W(r))4>^('r)  - E_^(r)  , (4) 

k k k 

where  W(r)  is  the  pseudopotential, 


W (r ) = Vc(r)  + J (E  - E ) | n«.><n£  | , (5) 

n,  1 k ,u 

and  the  pseudo-wave  function 

<J^<r)  = l a (k)  |k  + q>  (6) 

k q 1 

is  simply  a sum  of  plane  waves.  The  expansion  coefficients  a (k)  are 

q 

the  same  as  previously  used  in  the  expression  for  the  conduction 

eigenstate  [i.e.,  Equation  0]  . E^  is  the  eigenvalue  of  the  n,£  core 

c 

electron  state  in  the  crystal  environment.  Although  in  principle  t//  (r) 

k 

and  | n should  be  eigenstates  of  the  same  crystal  Hamiltonian  H,  in 
practice  |nfc>  is  well  approximated  by  an  eigenstate  of  the  free  ion. 

It  was  Phillips  and  Kleinman  (1959)  who  first  noted  the  pseudopotential 

c 

W(r)  is,  in  general,  weaker  than  the  true  crystalline  potential  V (r) . 

This  suggests  the  use  of  perturbation  theory  to  solve  for  the  eigen- 

states  <t>^(r)  and  the  eigenvalues  E of  the  pseudopotential  equation  (4). 
k k 

Note  that  the  eigenvalues  of  Equation  (4)  and  the  original  crystalline 
Schrodinger  equation  are  identical.  The  following  sections  will  out- 


line how  the  pseudopotential  equation  can  be  solved. 
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There  is  no  unique  solut ion  to  Equation  (4)  since  any  core-electron 

eigenstate  [nis  may  be  added  to  ^(r)  and  the  combination  will  still 

k 

♦ 

remain  a solution  of  the  pseudopotential  equation.  If  4'>lt)  could  be 

k 

obtained  exactly  (i.e.,  to  all  order  of  perturbation  theory),  this  non- 

-> 

uniqueness  would  not  matter  since  any  exact  solution  of  | (r)  would 

k 

•> 

yield  the  same  conduction  eigenstate  4'  (r) . This  is  evident  from  the 

k 

c ► ► 

relationship  between  4’,  (r)  and  d'  ^ C f ) shown  in  liquation  (7). 

k k 


U>,(r)  = (1  - P)<i>Jr) 
k k 


(7) 


Since  d (r)  will  be  obtained  to  first  order  in  perturbation  theory 

k 

it  is  necessary  to  further  specify  how  4*  (r)  can  be  made  unique. 

k 

Cohen  and  Heine  (1961)  have  suggested  that  an  appropriate  condition 

for  uniqueness  is  that  the  pseudopotential  equation  should  be  solved  for 

*■ 

the  'smoothest*  wave  function.  They  impose  the  criteria  that  4*  (r) 

k 

should  minimize  the  following  quantity, 


|V4>^(r)  | 2d3r/ 


I*  <r)|2d3i 


(») 


if  it  is  to  be  the  'smoothest*  wave  function.  This  condition  in  Equation 
(9)  can  be  shown  (Harrison,  i'H'O)  to  yield  the  following  equivalent 
criterion  for  the  smoothest  pseudo  wave  function. 


*nt  | W 1 ^ |4'„\  ^ J4V 

k k k k k k 


(9) 


Utilizing  Equ.it  ion  (9)  in  Equation  (-1)  leads  to  the  following  expression 
for  the  pseudopotential 


(l-P)v  |4vrvVr) 

W(r'  41  k (r ) = (1  - P)Vc4k(r)  + — - ~ 

k k - <*Jp|$>n 

k k k k 


(101 
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Squat  ion  (10)  is  the  nonlinear  optimized  pseudopotential  ti.e.,  Equation 

► 

(10) predicts  the  smoothest  0 (r)  form  of  the  pseudopotent ial ] . Said  other- 

k 

* * 

wise  the  optimized  form  of  W(r)  will  produce  a uniquely  defined  if  (r) 

k 

tliat  will  be  smoothest  m the  sense  of  equation  (H)  . 

Is  is  seen  that  the  optimized  pseudopotential  in  Equation  (10)  is  a 

► 

nonlinear  operator  since  it  depends  on  4*  (r) . However,  it  can  be  shown 

k 

(Harrison,  19t>6)  that  W(r)  may  be  linearized  without  affecting  the 

second  ordei  perturbation  theory  results  tor  t tie  eigenvalues  K . The 

k 

linear  optimized  form  of  W(r)  is, 


► I Cl  > 

W ( r ) if  (r)  - ().  - P)V^f  (r)  + — - V-_ Lk_  p (f  (1>> 
k k 1 - <k|p|k''  k 


(11) 


One  might  also  note  that  the  linearized  pseudopotent ial  in  Equation  (11) 

is  non-Hermitian  due  to  the  term  PV  . However, the  eigenvalues  are 

k 

- ► 

not  affected  by  the  non-hermit icity  of  W(r)  until  third  order  or  higher 
terms  are  reached  in  pertubation  theory  (Harrison,  196b). 

Perturbation  theory  predicts  that  the  eigenvalues  are  to  second 

k 

order , 


ii*k‘ 


t v k W 1 k ♦ 


k ♦ q|w!k'>'-k[wlk  + 


y ' K. * ip!  WJK.  ' ' 


g 


fr 

2m 


k‘ 


k + 


ql 2 


al 


(12) 


where  the  matrix  elements  of  W(r)  have  been  taken  between  plane  wave 
states,  and  the  sum  over  q has  tlie  term  q = 0 omitted.  With  a little 
algebra  and  noting  that  P‘  P,  and  Vv  = H - T,  it  can  be  shown  that  the 
diagonal  and  off-diagonal  matrix  elements  of  W are  (Harrison,  1966), 

(frk 


■ k | W | k ' 


-k | VC | kN  ♦ l 

n( 


2m  + ^k|V  |k'  - En« 


k|njj,^<njt|kN  (13) 
1 - sk|p|k- 


and 
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<k+q|w|k'*  = <k+q  | VC  | k^+  Y 


ft2k2 


-*■  | c I 

- E + <k | V | k> 


2m  nil 


y -e  5+<kivcik% 
‘ I 2n  rA  11 


n£i- 


<k+q  I nil'><nJl  | k>  + 

<k  | ni  "><nH  | k><k+q  | P | k>/  (l-<k|p|k>) 

(14) 


Equations  (13)  and  (14)  represent  the  unscreened  or  bare  pseudopotential 
matrix  elements.  The  matrix  elements  must  be  screened  before  they  are 
used  in  equation  (12).  The  screening  of  the  pseudopotential  is  dis- 
cussed in  later  sections  and  in  Appendix  D. 

Factorization  of  the  Pseudopotential 

Equations  (12),  (13)  and  (14)  can  be  used  to  evaluate  the  energy  of 

the  conduction  electrons  to  second  order  in  perturbation  theory,  but 

this  is  a difficult  task  because  of  the  complexity  of  evaluating  the 

crystalline  potential.  The  problem  of  computing  the  conduction-electron 

c 

energy  can  be  greatly  simplified  if  one  makes  the  assumption  that  V and, 
subsequently, W can  be  written  as  a linear  superposition  of  spherically 
symmetric  potentials  centered  upon  the  lattice  sites.  With  this  assump- 
tion, it  is  possible  to  factorize  the  piseudopotential  matrix  elements  in 
the  following  manner: 


k+q | W | kN  = S (q) <k+q | w | k>  , 


(15) 


where 


S (q)  = 1/N  \ exp  (-iq*r  , ) 

j 


and 


<k+q | w | k''  = 1/fi.  exp  [-i  (k+q)  *r  ] w(r)exp  (ik*r)  d r 


(16) 


(17) 
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The  quantity  S(q)  is  known  as  the  structure  factor,  r ^ is  the  position 

of  the  j’th  lattice  site,  il0  is  the  primitive  unit  cell  volume,  and 

w(r)  is  the  pseudopotential  associated  with  one  unit  cell. 


If  Equation  (15)  is  substituted  into  Equation  (12)  and  all  occupied 

■¥ 

k states  are  summed  over,  one  obtains  for  the  second  order  contribution 
to  the  conduction-electron  energy 

t ^ 

E = l | S (q) \2  F(q)  . (18) 

q 

F (q)  is  the  energy  wave  number  characteristic  and  is  explicitly 


F(|q|) 


2(1 


(270  3 


<k  | w ! k+q>«-k+qj  w|  k~*  djk 
^ (k3- | k+q| 2 ) 


•0  sc  + sc 
— n (q)  w (q) 


(19) 


The  singularity  in  the  denominator  of  the  first  term  of  Equation  (19) 
can  be  integrated  analytically  and,  therefore,  does  not  cause  any  great 
difficulty  when  |k|  = |k+q|. 


The  second  term  in  Equation  (19)  is  the  conduction  or  screening- 

electron  self-energy  which  has  been  subtracted  once  from  the  expression 

of  F(q)  because  the  first  term  has  counted  it  twice.  The  quantities 
SC  sc 

n (q)  and  w (q)  are  the  Fourier  transforms  of  the  screening-electron 
charge  density  and  Coulomb  potential,  respectively.  Equation  (18) 
represents  the  conduction-electron  energy  dependence  upon  the  crystal 
structure.  Hence,  if  the  crystal  is  distorted,  the  function  S(q)  changes 
and  affects  the  value  of  the  conduction-electron  energy  predicted  by 
Equation  (18).  Consequently,  Equation  (18)  is  of  critical  importance 
in  calculating  phonon  spectra  or  elastic  shear  constants  which  depend 
upon  a knowledge  of  the  distorted  lattice  energy.  Appendices  A and  B 
derive  the  explicit  relationship  between  F(q)  and  the  phonon  spectra  and 
elastic  shear  constants. 
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IV.  THE  CRYSTALLINE  POTENTIAL  AND  THE  PSEUDOPOTENTIAL 


The  derivation  of  the  HEP  pseudopotential  matrix  elements  was  out- 
lined in  the  precedinq  sections.  It  was  also  demonstrated  in  Appendices 
A and  B how  the  phonon  dispersion  curves  and  the  elastic  shear  constants 
of  a solid  metal  may  be  calculated  once  the  pseudopotential  matrix  ele- 
ments are  known.  Equations  (13)  and  (14)  indicate,  however,  that  there 
are  three  things  which  must  first  be  provided  before  a numerical  eval- 


uation of  the  HEP  pseudopotential  matrix  elements  can  be  made.  They  are: 

1)  the  core-electron  eigenfunctions  |n£>;  2)  the  core  electron  eigen- 

c 

values  E^^and  3)  the  crystalline  potential  V (r).  It  is  far  from  a 
trivial  task  to  obtain  any  one  of  these  three  quantities.  Therefore,  we 


shall  discuss  in  some  detail  the  analysis  for  obtaining  nH>,  E , V (r) 

n i 


and  the  approximations  that  have  been  employed. 


The  Core  Electron  Eigenstates 


Ideally, one  would  like  to  use  core-electron  wave  functions  and 
eigenvalues  that  have  been  calculated  self-consistently  for  the  crystal 
Hamiltonian.  However,  this  is  a formidable  if  not  impossible  task,  and 
one  must  resort  to  approximate  descriptions  for  the  core-electron  states 
in  the  crystal.  In  the  simple  metals  like  the  alkalis  where  there  is 
very  little  overlap  between  adjacent  atoms  in  the  crystal,  the  small 
core  approximation  may  be  utilized.  This  approximation  states  that  the 
core-electrons  in  the  crystal  are  described  by  the  wave  functions  for 
the  electrons  in  an  isolated  atom  or  ion.  A self-consistent  field 
calculation  was  performed  to  obtain  the  wave  functions  and  eigenvalues 
for  an  isolated  ion.  These  eigenstates  were  then  used  to  approximate 
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the  core-electron  states  in  the  crystal  for  the  purposes  of  constructing 
a first  principles  pseudopotential.  The  Herman-Ski 1 Iman  (1963)  atomic 
structure  program  was  used  to  (US  the  actual  numerical  computation  of 
the  free-ion  wave  functions.  The  Herman-Ski 1 Iman  (US)  program  performed 
a self-consistent  Hartree-Fock-Slater  determination  of  the  ionic  wave- 
functions.  One  could  also  do  a variational  determination  of  the  core 
states,  but  the  HS  self-consistent  field  approach  is  more  in  keeping 
with  the  first-principles  spirit  of  our  calculations. 

The  HS  program  also  calculates  the  free  ionic  eigenvalues  which  will 

be  denoted  as  E**\.  The  core-electron  eigenvalues  in  the  crystal,  K „ , 
tu  n x. 

will  have  different  values  because  of  interactions  with  the  crystal 
environment.  The  shift  in  the  core-electrons  energy  in  going  from  the 
free  ionic  to  the  crystalline  state  is  called  the  'core  shift'  (Lin  and 
Phillips,  1965).  As  will  be  shown  later,  it  is  important  to  do  an 
accurate  determination  of  file  core  shifts  when  doing  HFP  pseudopotential 
calculations.  The  accuracy  of  the  core  shift  can  significantly  affect 
the  calculated  phonon  dispersion  curves,  and  other  lattice  dynamical 
p roper t ies . 

The  crystalline  Potential 

The  crystalline  potential  was  constructed  as  a linear  sum  of 
spherically  symmetric  potentials.  Each  of  these  spherically  symmetric 
potentials  arises  from  the  charge  density  within  a single  primitive 
unit  cell.  The  potential  within  a primitive  unit  coll  'seen'  by  the 
conduction-electrons  consists  of:  1)  the  Coulombic  potential  of  the 
bate  i >n,  whore  the  ion  is  composed  of  a delta  function  nucleus  and  the 
localized  core-elect ron  charge  density,  2)  the  exchange  interaction 


between  the  core  and  conduction-electrons;  3)  the  orthogonalization 
hole  (Oil)  potential;  and  4)  the  self-consistent  screening  potential 
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1 


of  the  conduction-electron  gas. 


The  OH  and  screening  potentials  both 


arise  from  the  conduction-electron  charge  density  which  is  obtained 

c 

by  taking  the  absolute  square  of  i/i  (r)  and  summing  over  all  occupied 
k-states . 


c 

Contributions  2,  3 and  4 to  V are  discussed  in  the  next  two 
sections  and  in  Appendices  C and  D.  The  ion  potential,  contribution 
1,  as  previously  stated,  is  constructed  from  the  core-electron  eigen- 
states. The  eigenstates  yield  the  core  charge  density  and  then,  with 
the  aid  of  Poisson's  Equation,  one  obtains  the  core  potential. 


I 
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V.  APPLICATION  OF  THE  PSEUDOPOTENTIAL  TO  THE  CALCULATION 
OF  PROPERTIES  OF  THE  SOLID  ALKALIS 

The  accuracy  of  a given  pseudopotential  may  be  checked  by  using 
it  to  predict  properties  of  the  metal  that  can  be  verified  experi- 
mentally. The  elastic  shear  constants  and  the  phonon  spectra  are  two 
features  that  may  be  readily  predicted  with  the  aid  of  the  pseudo- 
potential and  compared  to  experiment.  The  phonon  dispersion  curves 
will  be  discussed  first  because  they  constitute  a more  rigorous  test 
of  the  correctness  of  the  energy  wave  number  characteristic.  The 
phonon  spectra  allows  a comparison  of  a set  of  curves  instead  of  a 
single  number  like  the  elastic  constants.  The  elastic  constants  depend 
on  F(q)  and  its  derivatives  only  at  the  reciprocal  lattice  vectors  G, 
whereas  the  phonon  frequencies  ui(q)  depend  on  the  values  of  F ( | G+q | ) . 
Hence  the  phonon  dispersion  curves  'sample'  a more  extended  range  of 
the  argument  of  F(q)  and  are  therefore  a more  sensitive  test  of  the 
correctness  of  F(q). 

Appendix  A gives  a detailed  mathematical  treatment  of  the  theory 
of  the  phonon  spectra  and  how  the  frequencies  u)(q)  are  calculated  once 
F(q)  has  been  computed.  In  this  section  we  present  a discussion  of  the 
sensitivity  of  the  phonon  spectra  to  certain  features  of  the  crystalline 
potential  used  in  the  construction  of  the  pseudopotential.  More  specif- 
ically, we  will  demonstrate  the  phonon  dependence  on  the  conduction- 
core  exchange  interaction,  the  core-electron  states,  the  conduction- 
electron  gas  screening,  and  the  shift  of  the  core-electron  eigenvalues, 


bL 


i.e.,  the  core  shifts. 
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Phonon  Spectra  a nd tt  i e Co  nduc  t i (ji  i-i'oi  e E xcl  lange  lnt  or  act  i on 

There  are  several  different  approximations  available  for  describing 
the  conduct ion- core  (C-C)  exchange  interaction.  Perhaps  the  best  known 
approximation  is  the  Slater  generalized  method  which,  in  the  case  of 
a nonuniform  charge  distribution,  becomes  the  so-called  local  density 
approximation  (LDA) 

v"  (r)  = -3  ( 3/rr ) 1 7 3 anl/,(r)  . (20) 

ex  c 

In  our  calculations,  n (r)  is  the  computed  core-electron  number  density. 
The  coefficient  a is  known  as  the  exchange  parameter  and  can  assume  a 
range  of  values.  In  our  analysis,  a was  chosen  to  have  one  of  the  three 
possible  theoretical  values  that  have  some  first-principles  justifi- 
cation. Slater  (1951)  originally  derived  Equation  (20)  for  the  case  of 
a free  electron  gas  for  which  the  value  of  a is  1.0.  Kohn  and  Sham  (1965), 
using  a variational  technique,  derived  an  expression  for  the  exchange 
potential  that  was  identical  to  Slater ' s, except  the  value  they  obtained 
for  a was  2/3.  More  recently,  I.indgren  and  Schwarz  (1972)  proposed 
using  values  of  a between  2/3  and  1.0  which  predict  total  energies  for 
a free  atom  that  are  equal  to  t tie  energy  predicted  by  an  exact  Hartree- 
f'ock  calculation. 

A second  approximation  for  the  C-C  exchange  potential  was  proposed 
by  Lindgren  (1971).  Hi'  suggested  .in  exchange  potential  given  by 

V1'  (r)  = -3  (3/tt)  1/3  [n'^lr)  - nl/3(r)l  , (21) 

ex  tv' 

where  the  total  electronic  density  n^  is  the  sum  of  the  conduction 

* -*  L 

[n  (r)|  and  core-electron  densities  In  (r)|.  Lindgren  derived  V 
v c ex 
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by  noting  that.  V(  is  more  properly  applied  to  the  total  charge  density, 

and  that,  if  one  wishes  to  study  the  interelectronic  conduction-core 

exchange, then  Equation  (21)  is  the  more  appropriate  form  for  describing 

this  interaction.  Note  that  the  exchange  potentials  in  Equations  (20) 

and  (21)  are  mite  similar  in  the  core  region  where  n >>  n , but  that 

c v 

V1  goes  to  zero  much  faster  in  the  region  where  n < n . Therefore 
ex  c v ' 

the  I.indgren  exchanqe  potential  predicts  a weaker  exchange  interaction 
at  large  r where  n , < nv  than  does  the  Slater,  Kohn-Sham,  or  hindgren- 
Schwarz  potentials. 

It  is  possible,  by  adjusting  the  value  of  a in  Equation  (20)  or  by 
introducing  a parameter  a in  Equation  (21),  to  fit  the  calculated  phonon 
spectra  to  the  experimental  data.  This  is  due  to  tire  fact  that  changing 
the  value  of  ot  in  a pseudopotential  calculation  produces  a shift  in  the 
magnitude  of  the  calculated  frequencies  without  changing  the  qualitative 
shape  of  the  curves.  However, the  temptation  to  paramaterize  it  has  been 
resisted  because  it  obscures  the  significance  of  other  contributions  to 
the  pseudopotential  (e.q.,  such  as  the  choice  of  dielectric  screening 
function).  Since  our  goal  of  a first  principles  pseudopotential  cal- 
culation was  to  try  to  elucidate  knowledge  about  fundamental  inter- 
actions in  the  solid,  all  calculations  in  this  thesis  utilized  one  of 
the  conduction-core  exchange  schemes  that  had  some  a priori  justifica- 
tion (e.q.,  Kohn-Sham,  Slater,  I.indgren,  and  Lindgren  and  Schwarz). 

To  illustrate  the  sensitivity  of  the  phonon  spectra  to  the  choice 
of  conduction-core  exchange  approximation,  we  show  in  Figure  1 two  phonon 
dispersion  curves  for  hi  that  were  calculated  using  identical  pseudo- 


potentials except  for  the  choice  of  C-C  exchange  potential.  The  solid 


curve  employed  VL 
ex 
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, and  the  dashed  curve  used  with  Kohn  and  Sham's 

ex 

value  for  a.  The  curves  are  chosen  only  to  illustrate  the  dramatic 

differences  produced  in  the  phonon  spectra  by  different  treatments  of 

the  conduction-core  exchange,  and  were  not  chosen  to  represent  the 

best  agreement  with  experiment.  The  dashed  curve  would  have  been  even 

a 

lower  if  Slater's  a was  used  for  Vgx.  Na,  K,  Rb  and  metals,  such  as 
the  hep  divalent  series,  show  similar  sensitivity  to  the  conduction- 
core  exchange. 


It  is  possible  to  make  some  general  observations  concerning  the 

conduction-core  exchange  approximations  of  Equations  (20)  and  (21) . 

Ot 

In  general  it  was  found  that  V of  Equation  (20)  predicted  phonon 

frequencies  for  Na  and  K that  had  the  longitudinal  acoustic  (LA)  and 

transverse  acoustic  (TA)  branches  in  the  (001)  direction  inverted  with 

respect  to  experiment.  The  inversion  of  the  LA  and  TA  modes  in  the 

(001)  direction  was  not  corrected  by  the  use  of  different  dielectric 

Ot 

screening  functions,  and  seems  to  be  a sympton  of  Vve  note  that  in 

the  other  crystal  directions,  (110),  (111),  where  the  LA  and  TA  modes 

ot 

are  not  so  nearly  degenerate,  V produces  good  agreement  with  experi- 
ment, with  the  LA  and  TA  branches  properly  ordered.  On  the  other  hand, 
the  Lindgren  exchange  approximation.  Equation  (21)  , does  predict  the  LA 

and  TA  branches  in  the  correct  order  for  the  (001)  direction  of  NA  and 

L d 

K.  Hence, for  Na  and  K,  the  Lindgren  potential  V ^ is  preferable  to  Vgx 
because  it  predicts  the  correct  ordering  of  the  phonon  modes  in  the 
(001)  direction,  as  well  as  yielding  better  qualitative  agreement. 


L Ot 

There  is  another  reason  why  V is  preferable  to  V when 

1 ex  ex 

culating  phonon  frequencies  for  Na  and  K.  Lindgren  exchange,  V 


cal- 

L 

ex 


appears  more  compatible  with  the  dielectric  screening  functions  which 

a 

go  beyond  the  llartree  approximation.  The  use  of  V in  conjunction 
with  the  Hartree  free  electron  gas  dielectric  function  (Sham  and  Ziman, 


1963)  predicts  phonon  frequencies  in  good  agreement  with  experiment 
(except  for  the  mode  inversion).  However,  when  the  interacting  elec- 
tron gas  dielectric  functions  are  used  with  V^x  or  , the  LA  phonon 

frequencies  are  lowered.  This  improves  agreement  with  experiment  if 

V1'  is  used,  but  further  aggravates  the  mode  inversion  of  the  (001)  LA 
ex 

and  TA  branches  which  occurs  with  V**  . The  dielectric  function  of 

ex 

Singwi,  et  al.  (1970)  (SSTL)  is  an  example  of  an  interacting  electron 
gas  screening  function  that  produces  this  sort  of  behavior. 


Although  definitely  produces  better  phonon  spectra  for  Na  and 

a L 

K than  does  V , the  same  cannot  be  said  for  Li.  V yields  phonon 
ex  ex 

frequencies  for  Li  that  are  much  too  high  in  comparison  to  experimental 

results,  regardless  of  the  dielectric  screening  function  that  was 

employed.  When  VU  is  used  with  the  Hartree  dielectric  function  one 
ex 

obtains  phonon  frequencies  for  Li  which  are  in  quite  good  agreement 
with  experiment  (seev  for  example,  Figure  6).  Furthermore,  the  cross 
over  of  the  La  and  Ta  branches  in  the  (001)  direction  of  Li  is  cor- 
rectly predicted.  In  the  calculations  presented  here,  the  quantity 
nv(r)  in  V^‘  was  approximated  by  l/i20  , the  uniform  conduction- 
electron  density.  Recent  calculations  by  Sun  (1978)  show  that  by 
approximating  nv(r)  in  \A  by  the  more  accurate  single  OPW  density 
[i.e.,  nv(r)  = z*A10  + n°*\r)],  the  phonon  curves  for  Li  were  sliqhtly 
lowered.  Though  this  was  still  not  enough  to  produce  good  agree- 
ment with  experiment , it  may  indicate  that  V1'  would  work  for  Li  if 

ex 


nv  was  known  sufficiently  accurately.  Sun  also  used  the  single  OPW 
result  for  nv  in  calculations  for  Na,  K,  and  Rb  and  found  that,  in  all 
cases,  the  more  rigorous  treatment  of  n^  in  improved  agreement 

with  experimental  phonon  curves. 

Phonon  Spectra  and  t he  Core- Electron  Eigenstates 

In  the  HFP  pseudopotential  formal  ism,  the  core-electrons  are  rep- 
resented by  either  an  analytic  or  numerically  computed  wave  function. 

This  sets  the  HFP  formalism  apart  from  model  pseudopotential  calculations 
which  do  not  explicitly  include  core-electron  states  in  the  construction 

of  the  pseudopotential.  The  Merman-Ski  1 lman  (1963)  atomic  st.ruc- 

HS 

ture  program  was  used  to  calculate  the  ionic  eigenvalues  (E^) > anc*  the 

radial  part  of  the  core-electron  wave  function  P „ (r)  for  the  ion  in  the 

nx. 

free  state.  The  MS  program  is  a self-consistent  field  calculation  based 
upon  a Hart ree-Fock-Slater  formalism. 

There  are  several  aspects  of  the  MS  program  that  deserve  some 
attention.  The  first  j>oint  is  that  the  core-electron  states  can  be  cal- 
culated with  different  approxi mat  ions  for  the  exchange  interaction  among 
the  core-electrons.  It  was  found  that  the  phonon  spectra  were  signifi- 
cantly affected  by  the  use  of  different  exchange  approximations  in  the 
HS  core  state  calculation.  In  Figure  2, there  are  two  phonon  dispersion 
curves  for  K that  were  calculated  with  identical  pseudopotentials  except 
that  the  solid  curve  used  core  states  calculated  with  Slaters  exchange 
[i.e.,  Equation  (20)  with  a equal  to  1]  and  the  dashed  curve  used  core 
states  calculated  with  Lindgren  and  Schwarz's  value  for  the  a-parameter 
[Equation  (20)  with  a equal  to  0.722],  As  can  be  seen, there  is  a very 
significant  difference  in  the  two  curves.  It  is  tempting  to  conclude 


*1* 


the  construction  of  the  core-electron  states.  Both  curves  were  computed  from 
identical  pseudopotentials  except  the  solid  curve  resulted  from  core  states 
calculated  with  Slater’s  exchange  approximation  [a  = 1.0  in  Equation  (20)] 
and  the  dashed  curve  from  Lindgren  and  Schwarr/  approximation  [a  = 0.722  in 
Equation  (20)].  Experimental  data  points  for  K are  from  Cowley  et  al.  (1966) 
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on  the  basis  of  Figure  2 that  the  Lindqren  and  Schwarz  exchange  approxi- 
mation provides  a better  description  of  the  core-electron  states,  but 
as  we  shall  see, this  is  not  necessarily  so.  Further  study  of  the  effect 
of  the  exchange  approximation  upon  the  core  states  reveals  that  Slater's 
approximation  yields  core  state  eigenvalues  that  are  in  better  agree- 
ment with  experiment  than  eigenvalues  predicted  by  any  of  the  other 
exchange  approximations  previously  discussed.  Therefore,  using  as  a cri- 
terion the  difference  between  the  experimental  and  theoretical  eigen- 
values, one  finds  that  Slater's  approximation  for  the  exchange  inter- 
action is  the  most  appropriate  choice  for  the  core  wave  function  calcu- 
lation. Thus,  the  apparently  excellent  agreement  in  Figure  2 predicted 
by  the  Lindqren  and  Schwarz  core  states  may  bo  partly  fortuitous. 

Although  the  use  of  Slater's  exchange  in  the  construction  of  the  HS 
core  states  gives  best  agreement  between  experimental  and  theoretical 
core  eigenvalues,  differences  that  remain  are  still  significant.  Day 
et  al.  (1976)  showed  for  Li  and  Be  that  the  differences  between  the 
theoretical  and  experimental  core  eigenvalues  could  measurably  affect 
the  resulting  theoretical  phonon  dispersion  curves.  Figure  3 illus- 
trates, for  Li,  the  differences  in  the  phonon  spectra  produced  when 
theoretical  and  experimental  core  eigenvalues  are  used.  In  principle 
the  experimental  ionic  eigenvalues  (with  appropriate  core  shifts' 
should  always  be  used  in  HFP  pseudopotential  calculations,  out  for  multi- 
orbital cores  like  Na  and  K,  there  is  little  reliable  spectroscopic  data 
available  for  the  ion.  Theoretical  term  values  were  employed  in  our 
calculations  due  to  the  unavailability  of  experimental  data  for  the 
inner  shell  ionic  energy  levels.  Table  1 lists  experimental  and  theoret- 


*4 


ical  eigenvalues  for  the  Li  ion  and  the  Na  and  K atoms  [Slater  (1955)  , 
Bearden  and  Burr  (1967),  CRC  Handbook  (1972)]. 
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Table  1.  Experimental  versus  theoretical  eigen- 


values for  the  Li+  ion  and  the  Na  and 
K atoms . 


Ene 

Experiment 

Theory 

( Rydbergs ) 

( Rydbergs) 

. + 

Llls 

- 5.56 

- 5 . 38 

Na, 

Is 

- 79.4 

- 78.0 

2s 

- 5.2 

- 4.7 

2p 

- 2.8 

- 2.67 

3s 

- 0.38 

- 0.37 

- 

Is 

-266.2 

-262.1 

2s 

- 28.2 

- 27.1 

2p 

- 22.2 

- 22.0 

3s 

- 3.0 

- 2.95 

3p 

- 1.81 

- 1.73 

4s 

- 0.32 

- 0.308 

The  Conduction-Electron  Potential  and  the  Phonon  Spectra 

The  conduction-electron  charge  density  is,  for  convenience, 

separated  into  two  different  parts,  the  orthogonalization  hole  density 

OH  ♦ sc  * OH  ♦ 

n (r)  and  the  screening  charge  density  n (r).  The  density  n (r) 

-► 

is  independent  of  the  perturbing  potential  W(r)  and  is  the  result  of 
orthogonal izing  the  plane-waves  to  the  core  wave  functions.  The 
screening  charge  represents  the  linear  response  of  the  conduction- 
electron  gas  to  W(r).  These  charge  densities  each  give  rise  to  their 
own  respective  potentials  which  are  important  contributions  to  the 
crystalline  potential.  They  also  interact  with  the  core-electron 

charge  density  which  produces  a shift  in  the  eigenvalues  of  the  ^ore- 

OH  o pw 

states.  The  core  shift  due  to  n has  been  given  the  special  name  v f . 

Cutler  et  al.  (1975)  and  Sun  et  al.  (1976)  have  shown  in  calculations 

of  phonon  spectra  that  it  is  important  to  do  an  'exact'  treatment  of 
OH  . , , opw  , 

n (r)  and  v , instead  of  the  approximate  treatments  proposed  and  used 
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by  Harrison  (1966).  A detailed  derivation  and  discussion  of  the 

exact  and  approximate  OH  treatments  is  presented  in  Appendix  C. 

Moriarty  (1970)  has  used  an  exact  treatment  of  the  OH  for  the  noble 

metals  Cu,  Ag,  and  Au  but  the  calculations  in  this  thesis  are  the 

OH 

first  to  use  an  exact  n for  the  alkali  metals. 

The  effect  of  the  exact  treatments  of  n°H  and  v0t>W  is  seen  in 
Figure  4 where  the  phonon  spectra  for  Li  calculated  with  the  exact 
OH  (solid  curves)  is  compared  with  a similar  calculation,  using,  how- 
ever, an  approximate  OH  (dashed  curves) . The  effect  of  the  exact 
OH 

n is  to  significantly  lower  (i.e.,  by  5 to  10%  at  the  Brillouin  zone 
boundaries)  the  phonon  dispersion  frequencies.  Similar  behavior  is 
noted  for  Na  and  K. 

SC 

A derivation  of  the  screening  potential  w is  given  in  Appendix 
D,  where  it  is  pointed  out  that  the  screening  of  the  conduction-electron 
gas  can  be  described  by  a scalar  or  static  dielectric  response  function 

6*(q) • 

e*(q)  = 1 + (l-E(q) ) (l-G(q) ) . (22) 


The  function  e(q)  is  the  free-electron  or  Hartree  dielectric  function 
[see  Equation  (a  10))  and  G(q)  is  a function  which  accounts  for  the 
effects  of  electronic  exchange  and  correlation  between  the  conduction 
electrons.  For  a free-electron  gas,  G(q)  is  zero  and  e*(q)  becomes 
equal  to  £(q).  Although  there  is  no  exact  solution  for  G(q)  for  an 
interacting  system, there  have  been  many  attempts  at  calculating 
approximate  forms  for  this  exchange-correlation  function  [Toigo  and 
Woodruff  (1971),  Geldart  and  Taylor  (1970),  Hedin  and  Lundqvist  (1971)  J. 
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For  a number  of  reasons  the  G(q)  of  Singwi  et  al.  (SSTL)  (1970)  was 
employed  in  most  of  our  calculations.  Firstly,  the  SSTL  evaluation 
of  G(q)  can  be  represented  by  tin?  simple  analytic  expression  of 
Equation  (23) 


G(q) 


A 


-Bq2/k 

1-e 


? 


(23) 


winch  is  convenient  for  use  in  numerical  calculations.  Secondly  the 
dielectric  function  must  satisfy  certain  relations  involving  the  pair 
correlation  function  g(r)  for  r *■  0 and  the  isothermal  compressibility 
(the  inverse  of  the  bulk  modulus)  SSTL  (1970).  The  SSTL  dielectric 
function  makes  reasonable  predictions  for  g(r)  at  small  r and  closely 
satisfies  the  so-called  compressibility  sum  rule.  The  limiting  values 
of  G (q)  go  from  zero  at  small  q to  - 1 at  large  q,  i.e.,  q > 2k  . 


Figure  5 shows  two  sets  of  phonon  spectra  for  Na  that  were  cal- 
culated with  identical  pseudopotentials  except  that  the  dashed  curve 
employed  the  SSTL  G(q)  and  the  solid  curve  was  calculated  with  the 
Hartree  t'(q),i.e.,  G(q)  = 0.  Significant  differences  are  observed  in 
the  two  sets  of  phonon  spectra,  though  not  as  significant  as  those  that 
are  produced  by  different  treatments  of  the  conduction-core  exchange 
interaction  or  the  core-electron  states.  The  effect  of  non-zero  G(q) 
dielectrics  is  to  reduce  the  ]J\  phonon  frequencies  with  respect  to 
those  predicted  by  the  Hartree  dielectric  function.  As  one  would 
expect,  the  TA  frequencies  are  not  as  strongly  influenced  by  the 
dielectric  screening  function.  Our  calculations  indicate  that  the 
specific  form  of  G(q)  is  not  critically  important;  however,  it  is 
important  to  include  some  form  of  the  exchange-correlation  screening 


in  the  calculations. 
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VI.  PHONON  SPECTRA  AND  ELASTIC  SHEAR  CONSTANTS 
OF  THE  SOLID  ALKALI  METALS 

Results  for  la,  Na  and_K 

In  the  previous  sections,  it  was  shown  how  an  HFP  pseudopotential 
is  formulated  and  how  the  crystalline  potential  can  be  constructed. 

In  addition,  the  dependence  of  the  phonon  spectra  on  various  contri- 
butions to  the  crystalline  potential  was  demonstrated.  More  specif- 
ically, we  considered  the  dependence  of  the  crystalline  potential  on 
the  conduct  ion-core  exchange  interaction,  ami  the  dielectric  screening 
of  the  conduction  electron  gas.  Pseudopotentials  were  calculated  for 
Li,  Na , and  K using  various  combinations  of  the  conduction-core  ex- 
change potential  (i.e.,  V*'  or,  V^x  with  one  of  the  possible  a priori 
values  of  ct)  and  the  dielectric  screening  function  (i.e.,  the  Hartree 
or  SSTL  dielectric  functions).  The  validity  of  each  pseudopotent ial 
was  tested  by  calculating  the  phonon  spectra,  elastic  shear  constants 
and  the  ion-ion  pair  potentials.  In  this  section  the  results  of  these 
calculations  will  be  presented. 

All  pseudopotentials  that  were  tested  included  an  exact  treatment 
of  the  oil  density.  The  core-electron  states  were  calculated  by  the  HS 
atomic  structure  program  which  has  been  discussed  previously  in  Sections 
IV  and  V.  For  Na  and  K,  it  was  found  that  a combination  of  Lindgren's  ex- 
change i>otential  and  the  SSTL  dielectric  screening  function  in  the 
pseudopotential  gave  the  best  overall  agreement  with  experimental  phonon 
spectra  and  elastic  shear  constants.  On  the  other  hand,  the  best  results 
for  Li  were  obtained  with  the  Kohn-Sham  description  of  the  conduction 
core  exchange  and  the  Hartree  dielectric  function  for  the  conduction 


electron  gas.  The  phonon  spectra  predicted  by  these  pseudopotentials 
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for  Li,  Na  and  K are  shown  in  Figures  6,  7 and  8,  respectively . Table 
2 lists  the  elastic  shear  constants  for  Li,  Na  and  K predicted  by  these 
same  pseudopotentials.  The  details  of  the  elastic  shear  constant  calcu- 
lation for  these  BCC  materials  are  given  in  Appendix  B.  For  the  latter, 
theoretical  results  are  generally  within  20%  or  better  of  the  experi- 
mental data. 

These  results  can  be  considered  as  quite  good  when  one  remembers 
the  fundamental  nature  of  the  a priori  pseudopotential  calculations. 

It  is  worth  noting  that  we  have  also  studied  one  parameter  schemes  which 
produce  pseudopotentials  that  yield  nearly  perfect  agreement  with  experi- 
mental phonon  spectra  and  elastic  shear  constants.  However,  the  use  of 
the  a priori  approach  is  capable  of  producing  insight  into  some  of  the 
fundamental  interactions  in  the  solid  which  the  model  or  paramaterized 
pseudopotential  cannot,  in  general,  do.  It  is  encouraging  to  note  that 
the  best  results  for  Na  and  K were  obtained  with  consistent  treatments 
of  the  conduction  core  exchange  interaction  and  conduction  gas  screening. 
In  these  calculations  for  Na,  and  K,  the  conduction  electron  density  nv<r) 
in  Equation  (21)  for  Lindgrente  exchange  was  approximated  by  the  average 
plane  wave  density  1/S10 . Sun  (1978)  has  shown  that  by  approximating 
nv(r)  by  the  1-OPW  density,  the  phonon  spectra  for  Na  and  K are  brought 
into  almost  perfect  agreement  with  experiment.  He  also  obtained  similar 
excellent  results  for  the  heavier  alkali  metal  Rb.  This  reinforces  the 
evidence  that  Lindgren's  C-C  exchange  approximation  is  a good  description 
for  the  alkalis.  It  is  disappointing  but  not  really  surprising  that 
the  pseudopotential  for  Li  required  a different  treatment  of  these  same 


A 


interactions.  Ham  (1972),  in  his  classic  calculations  of  the  band 
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Table  2. 

Elastic  shear 

constants  for  the 

simple 

> alkali! 

i Li,  Na,  and  K. 

Units 

of  1012 

dyne/cm‘ . 

C'l  1-C1  2 

Li 

exp 

0.108 

0.0230 

th 

0.126 

0.0187 

Na 

exp 

0.0577 

X 

O 

th 

o 

o 

0.0137 

exp 

X 

fN 

o 

d 

0.0075 

th 

0.0355 

0.0080 

structure  for  the  alkali  metals  found  that  Li  has  the  greatest  relative 
deviation  from  free  electron  like  behavior.  When  Lindgren  exchange  was 
used  for  Li,  it  was  found  that  the  phonon  frequencies  were  raised  much 
too  high  regardless  of  what  dielectric  function  was  used.  Although  we 
have  obtained  pseudopotentials  for  all  the  simple  alkalies  which  yield 
reliable  results  for  the  solid  state  and  Sun  has  shown  how  refinements 
may  be  made,  the  anomalous  behavior  of  Li  suggests  that  as  of  yet  there 
exists  no  single  unified  description  for  constructing  the  pseudopotential 
of  simple  metals. 
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VII.  CALCULATION  OF  THE  LIQUID  METAL  PROPERTIES  OF  THE  ALKALIS 

It  has  been  shown  in  the  preceding  sections  that  a first  principles 
non-local  pseudopotential  can  reliably  predict  properties  of  simple 
metals,  such  as  the  alkalis,  in  the  solid  state.  This  was  demonstrated 
explicitly  for  the  elastic  shear  constants  and  phonon  spectra  of  Li, 

Na,  and  K.  We  will  now  show  that  it  is  possible  to  predict  properties 
of  the  liquid  alkali  metals  using  a first  principles  pseudopotential 
with  appropriate  modifications  for  the  liquid  state. 

Liquid  Metal  Structure 

Before  determining  the  structure  of  liquid  Li,  Na,  and  K using 
pair  potentials  that  are  obtained  from  a pseudopotential  calculation, 
a brief  discussion  of  what  is  meant  by  'liquid  structure’  will  be 
presented . 

The  concept  of  structure  obviously  has  different  implications  in 
liquid  and  solid  systems.  The  particles  in  a solid  have  definite 
lattice  positions,  but  in  a liquid  the  particles  are  not  locked  into 
definite  locations  and  are  much  freer  to  move  about  than  in  a solid. 

The  motion  and  positions  of  the  particles  in  the  liquid  will  be 
correlated  because  of  the  interactions  between  the  particles.  Hence 
when  one  says  a liquid  has  structure,  one  implies  that  there  is 
correlated  behavior  of  the  particles  in  the  liquid  and  not  that  there 
is  a rigid  lattice  structure. 

Perhaps  the  most  direct  way  to  recognize  that  there  is  structure 


in  a liquid  system  is  to  observe  the  results  of  an  X-ray  or  neutron 
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diffraction  experiment  performed  upon  a real  liquid,  and  contrast 
this  with  the  ideal  liquid  and  ideal  crystal  cases.  In  any  diffrac- 
tion experiment,  the  object  is  to  obtain  the  diffracted  intensity 
1(0)  of  X-rays  or  neutrons  as  a function  of  the  scattering  angle  0. 

For  an  ideal  crystal  with  perfect  ordering, the  diffracted  intensity 
I is  zero  except  at  a finite  number  of  angles  where  sharp  spikes 
occur  in  the  intensity  of  the  diffracted  X-rays  or  neutrons  (see 
Figure  9) . This  is  the  well  known  Ewald  sphere  scattering  which  is 
discussed  in  any  elementary  solid  state  physics  text  [e.g.,  Ziman  (1972)]. 
Consider  now  an  ideal  liquid  which  has  perfect  disorder.  In  the  ideal 
liquid, there  is  no  correlation  between  the  positions  of  the  particles 
and  it  can  therefore  be  considered  to  be  a homogeneous  continuum.  The 
diffracted  intensity  will  show  no  angular  dependence  because  all 
angles  of  scattering  are  equally  likely  for  this  ideal  limit.  The 
ideal  liquid  would  display  a perfectly  uniform  diffracted  intensity, 
as  is  illustrated  in  Figure  9.  On  the  other  hand,  1(0)  for  a real 
liquid  does  show  a definite  angular  dependence  that  lies  somewhere 
between  the  discrete  angle  scattering  of  the  ideal  crystal  and  the 
diffuse  scattering  of  the  ideal  liquid  as  Figure  15  illustrates. 

One  can  obtain  from  experimental  diffraction  experiments  a 
function  S (q)  known  as  the  interference  function  or  the  structure 
factor.  S (q)  is  proportional  to  the  square  of  the  diffracted  in- 
tensity and  q is  the  wavenumber  of  the  diffracted  radiation  which  is 
related  to  the  scattering  angle  (March,  1966).  S (q)  can  also  be 
calculated  theoretically.  The  formal  definition  of  S (q)  is  given 


by 
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S (q ) 


t t 


Figure  9.  Structure  factor,  S(q),  curves  for  ideal  systems. 

Solid  curve  - liquid  system,  dashed  curve  - ideal 
crystal.  S(q)  is  proportional  to  the  square  of 
the  diffracted  intensity. 
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Figure  10.  Model  of  a two-dimensional  liquid  with  N = 5 
particles  per  cell  and  periodic  boundary 
conditions.  Each  particle  is  allowed  to 
interact  with  its  N-l  nearest  neighbors. 
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S(q) 


1/N 


l exp  liq* (r . -r  . ) ] 

i.  j 1 -1 


(26) 


where  the  i,j  sum  is  over  the  coordinates  of  the  N particles  and  the 
bar  denotes  that  an  ensemble  or  configurational  average  is  to  be 
taken  (Faber,  1972). 


Another  manner  in  which  the  liquid  structure  can  be  specified  is 
through  a hierarchy  of  correlation  functions  (Croxton,  1974)  . In 
general  one  may  define  2,3,...,  N-body  correlation  functions  which 
completely  specify  the  liquid  structure  if  all  are  known.  The  two- 
body  correlation  function  g(r)  is  also  known  as  the  pair  correlation 
or  radial  distribution  function.  The  function  g(r)  represents  the 
probability  density  that  a pair  of  particles  in  the  liquid  are  a 
distance  |r|  apart.  If  one  particle  is  at  the  origin, then  the  prob- 
ability that  a second  particle  is  in  d'r  about  r is  n0g(r)d3r,  where 
n„  is  the  average  number  density.  The  pair  correlation  function  is  of 
particular  interest  because  it  and  S (q)  form  a Fourier  transform 
pair : 


S (q) -1  = n„  | (g(r)-l)  exp 


To  be  more  precise,  S (q) -1  and  g(r)-l 
pair . 


( iq • r ) d 3 r 

constitute  the  Fourier 


(27) 

transform 


Closely  related  to  g(r)  is  the  cumulative  distribution  function 
G(r).  This  function  represents  the  total  number  of  neighbors  within 
a distance  r of  a particle  located  at  the  origin  r = 0.  In  terms  of 
the  pair  correlation  function,  G(r)  is 


G(r)  = 4tt  n0 


g (r ' ) r ' 2dr ' 


0 


(28) 


42 


G(r)  and  g(r)  specify  how  pairs  of  particles  are  distributed  in  the 
liquid,  and  in  general  they  are  related  to  three-body  and  higher 
correlation  functions;  these  will  not  be  discussed  in  this  work. 

Having  defined  the  meaning  of  'liquid  structure,'  one  is  now 
confronted  with  the  task  of  determining  that  structure.  The  inter- 
actions between  the  particles  of  a liquid  are  responsible  for  its 
structure.  If  the  particles  in  the  liquid  are  assumed  to  interact 
via  a spherically  symmetric  two-body  pair  potential  U(r),  then  the 
energy  of  the  j'th  particle  interacting  with  its  neighbors  in  an  N 
particle  system  is 
N 

E.  = I U(|r.-r.|)  . (29) 

1 • u ■ 1 3 

i/j 

The  effective  ion-ion  pair  potential  U (r)  can  be  calculated  from 
pseudopotential  theory  and  will  be  discussed  in  detail  later. 

Methods  for  Calculating  the  Liquid  Structure 

There  are  a number  of  theories  which  relate  the  potential  U(r) 
to  the  liquid  structure.  Perhaps  the  best  known  of  those  which  have 
had  some  success  in  the  past  are  the  Percus-Yevick  (1958)  and  the 
Born-Green  (1946)  theories.  However,  there  are  much  more  direct  and 
reliable  approaches  to  the  calculation  of  the  liquid  structure. 

These  are  the  Monte  Carlo  (MC)  and  Molecular  Dynamics  (MD')  computer 
simulations.  The  MC  and  MD  methods  use  the  computer  to  simulate  the 
actual  behavior  in  the  liquid  system.  In  the  MD  method,  one  solves 
the  equations  of  motion  for  N particles  in  a box  with  periodic 
boundary  conditions.  After  some  assumptions  have  been  made  for  the 


particle's  initial  positions  and  velocities  their  subsequent  trajectories 


can  be  computed  if  the  interaction  potential  t)(i)  is  known.  In  the  Ml) 


method,  the  computer  records  the  particles  coordinates  as  time  progr esses 


and  computes  a time  average  of  whatever  property  of  the  system  is  of 


interest  . On  the  other  hand,  the  Monte  Carlo  method  computes  an  en- 


semble average  value  of  f(x),  wher  e f (x)  can  be  any  static  property  of 


the  liquid  such  as  its  structure.  in  the  Mr-  calculat ion, the  computer 


generates  random  confiqurat  ions  of  the  N particles  and  calculates  the 


relative  probability  of  each  configuration's  occurrence.  The  ensemble 


average  value  of  fix)  is  formed  by  calculating  i t s conf igurat ion 


weighted  average.  The  next  tew  sect  ions  will  discuss  the  details  ot 


this  process. 


The  Monte  Carlo  Calculation  and  the  Markov  Chain 


The  Monte  Carlo  technique  can  be  used  t o esl imato  the  ensemble 


or  configurational  averages.  II  s " (>,»•  ,»  •••  > r ) is  tin*  N- 

1 ' N 


body  collective  coordinate  for  a system  ot  N interacting  particles 


confined  to  a volume  Si,  and  f (x)  is  any  non-singular  reasonable  func 


t ion  of  x.thon  the  pet it-canonical  ensemble  average  value  of  f txl  is 


f (x)  exp  l-|l  U (x)  lil  x 
N 


* ' N 

exp  l -|i  U ( x)  1 d x 
N 


where  t)  (x!  is  t hi'  interaction  potential  between  the  N particles,  and 

N 


(i  is  the  usual  Boltzmann  factor  1/KT.  As  is  well  known,  the  IN 


Kh 


dimensional  integrals  in  liquation  t '0)  have  analytic  evaluations  for 
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only  a very  small  class  of  U (x)  and  f(x) . Therefore,  it  is  in  general 

N 

necessary  to  have  a method  of  approximating  Equation  (30)  for  arbi- 
trary U and  f.  The  Monte  Carlo  method  calculates  the  ensemble  average 
value  of  f (x)  by  using  the  following  expression 


N 

c , 

f - 1/N  y f(x.)  . (31) 

c • , l 


The  sum  in  Equation  (31)  is  over  possible  configurations  x^,  that  the  co- 
ordinates of  the  N particles  may  assume,  and  N is  the  total  number  of 
configurations  included  in  the  sum.  The  prime  on  the  sum  denotes  that 
the  configurations  x^  included  in  the  sum  were  chosen  in  a particular 
way  to  constitute  a so-called  'Markov  Chain.'  Choosing  the  set  of 

configurations  (x.)1  so  that  they  constitute  a Markov  Chain  ensures 
l i=l 

that  each  x^  enters  the  sum  in  Equation  (31)  with  its  relative  prob- 
ability of  occurrence.  Therefore  Equation  (31)  represents  a simple 
weighted  average.  The  above  discussion,  although  brief,  sums  up  the 
essence  of  the  MC  technique.  A more  detailed  discussion  is  given  in 
Metropolis  (1953)  or  Wood  (1968) . 


The  liquid  structure  of  the  alkali  metals  was  determined  using 
the  Monte  Carlo  method.  Since  this  method  is  basically  a computer 
simulation  technique,  a model  of  the  liquid  had  first  to  be  established. 
The  liquid  metal  system  was  represented  as  a cube  of  side  L containing 
N particles.  The  values  of  N and  L were  chosen  so  that  N/L3  equalled 
the  average  number  density  of  the  liquid  system  being  studied.  The 
number  of  particles  in  the  cube  was  chosen  arbitrarily  as  either 


125  or  216.  Surface  effects  for  small  systems  of  this  size  could 
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appreciably  affect  the  structure.  To  limit  the  effect  of  the  surface 
upon  the  structure  and  to  simulate  an  infinite  system,  periodic  boundary 
conditions  were  employed.  Hence, if  a particle  in  the  box  moved  out 
through  one  face  of  the  cube, an  identical  particle  would  enter  the 
opposite  face.  The  particles  in  the  cube  interacted  with  each  other 
and  with  particles  in  adjoining  cubes  (adjoining  cubes  due  to  periodic 
boundary  conditions)  via  an  effective  ion-ion  pair  potential  U (r)  . As 
previously  mentioned,  U(r)  was  calculated  using  first  principles  pseudo- 
potential theory.  A two-dimensional  analog  of  this  three-dimensional 
system  is  illustrated  in  Figure  10  with  only  five  particles  in  each 
cell  for  the  sake  of  simplicity. 

The  generation  of  the  configurations  constituting  the  Markov  Chain 
was  accomplished  by  randomly  moving  the  particles  in  the  cube  one  at  a 
time.  The  move  of  a particle  was  accomplished  with  the  aid  of  a random 
number  generator.  Three  random  numbers  were  used  to  change  the  x,y, 
and  z coordinates  of  a particle  every  time  it  was  moved . Each  time  a 
particle  was  moved, a test  was  performed  to  determine  whethei  the  new 
configuration  would  be  retained  as  a member  of  the  Mat kov  Chain  or 
whether  it  would  be  rejected  and  the  original  configuration  would  be 
recounted  as  a member. 

The  exact  method  for  testing  whether  a configuration  will  be 
accepted  or  rejected  as  a Markov  Chain  member  proceeds  as  follows. 

A particle  is  moved  from  its  initial  position  to  a trial  position. 

The  energies  K1  and  E*  are  calculated  for  the  particle  in  its  initial 
and  trial  positions,  respectively.  These  energies  are  then  used  to 
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calculate  the  Boltzmann  probability  factors  P*  and  Pfc  corresponding 


to  the  initial  and  trial  positions,  where  P is  given  by 


P = exp  (-E/k  t)  • (32) 

The  energies  E1  and  E*"  are  computed  by  using  Equation  (29)  to  sum  over 
the  N-l  nearest  neighbor  pair  interactions.  The  inclusion  of  the  N-l 
nearest  neighbors  in  the  energy  sum  is  known  as  the  'minimum  image 
convention'  and  is  illustrated  for  the  two-dimensional  system  of 
Figure  10  (Wood  and  Parker,  1957).  When  P*"  1 P1,  the  trial  configura- 
tion is  automatically  accepted  as  a member  of  the  Markov  Chain  and 
it  becomes  the  initial  configuration  the  next  time  a particle  is 
moved.  However,  if  P1  > p'",  then  the  ratio  Pt/P1  is  formed  and  a random 
number  0 S r < l is  also  generated.  If  pVp1  - R,then  the  trial  con- 
figuration is  accepted  as  a Markov  Chain  member  and  again  becomes  the 
initial  configuration  for  the  next  move.  If  Pt/P1  < R,  then  the  trial 
configuration  is  rejected,  the  initial  configuration  is  recounted  as 
a Markov  Chain  member,  and  the  next  move  uses  the  retained  configura- 
tion as  its  starting  point. 

The  Markov  Chain  is  composed  of  many  configurations,  each  of 
which  is  generated  by  a random  particle  move  that  satisfied  the  above 
criteria.  The  maximum  allowed  size  of  the  random  moves  can  affect 
how  efficiently  the  Markov  Chain  is  generated.  If  the  moves  are  too 
big,  there  is  a relatively  high  probability  that  a large  percentage 
of  configurations  will  be  rejected  because  particles  will  be  moved 
too  close  together.  On  the  other  hand,  if  the  moves  are  too  small, 
then  each  configuration  will  differ  insignificantly  from  the  last  one 
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and  it  will  take  an  excessive  number  of  configurations  to  produce  an 

average  value  of  f(x)  in  lljuation  (31).  In  practice. we  choose  the 

l /, 

maximum  allowed  displacement  of  the  random  move  to  be  0.1(l/nf) 

This  choice  resulted  in  approximately  50%  of  the  generated  configura- 
tions being  included  in  the  Markov  Chain. 

The  Monte  Carlo  technique,  as  described  above,  was  used  to  com- 
pute the  ensemble  average  value  of  the  cumulative  distribution  func- 
tion G(r)  for  a discrete  set  of  r values.  This  function  was  evaluated 
for  each  configuration  in  the  Markov  Chain  as 

N 

G(r)  = l 0 (r- 1 r . -r  | ) , (33) 

i=l  1 3 

where  0 is  the  step  function,  and  the  sum  over  i counts  the  pair 
| “►  ”►  | 

distances  |r^-r  | between  the  j'th  particle  and  its  N-l  neighbors. 

The  prime  on  the  sum  denotes  that  only  the  N-l  nearest  neighbor 
pairs  were  included  in  the  evaluation  of  G(r),  and  r is  the  co- 
ordinate of  the  particle  just  moved  to  form  the  configuration  under 
consideration.  The  above  procedure  was  used  by  Brush  et  al.  (1966) 
in  the  evaluation  of  the  structure  of  a one-component  plasma.  As 
they  pointed  out, it  is  redundant  in  the  long  run  to  use  pair  distances 
in  Equation  (33)  which  do  not  involve  the  particle  whose  move  was 
just  accomplished. 

The  inclusion  of  only  the  N-l  nearest  neighbor  pairs  in  Equa- 
tion (33)  has  the  consequence  that  the  cumulative  distribution  function 
will  only  be  defined  on  the  interval  0 5 r - L/2.  One  could  go  beyond 


the  minimum  image  convention  and  evaluate  G(r)  for  distances  greater 
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than  L/2,  but  the  periodic  nature  of  the  system  and  its  finite  size 
would  make  the  accuracy  of  G(r)  at  these  distances  questionable.  To 
extend  G(r)  to  larger  r by  using  bigger  systems  (i.e., larger  N and  L) 
requires  an  inordinate  amount  of  computer  time  and  rapidly  becomes 
impractical . 
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VIII.  THE  PAIR  POTENTIAL 


The  pair  potential  U(r)  representing  the  effective  interaction 
between  the  ions  in  a metal  can  be  written  as  (Harrison,  1966) 


U(r) 


/r  + 


2S20/  (2tt) 


where  F(q)  is  the  energy  wave  number  characteristic  [Equation  (19)] 
and  2*  is  the  effective  valence  of  the  ions  [Equation  (C.8)].  Z* 
is  found  by  adding  up  the  nuclear,  core  electron,  and  orthogonaliza- 
tion  hole  charges.  The  first  term  in  Equation  (34)  represents  the 
direct  Coulomb  interaction  between  the  ions.  The  second  term  is  an 
indirect  interaction  representing  the  influence  of  the  conduction 
electron  gas.  The  indirect  interaction  is  obtained  from  the  second 
order  contribution  to  the  conduction  electron  energy  E which  is 
given  in  Equation  (18)  . E contains  the  absolute  square  of  the  struc- 
ture factor  S (q) . This  means  E will  include  an  explicit  summation 
over  pairs  of  atoms,  and  therefor^  each  term  in  the  pair  sum  in  the 
second-order-conduction  electron  energy  E can  be  identified  as  an 
effective  ion-ion  potential  interaction  contributing  to  the  total 
energy. 

The  F (q)  in  Equation  (34)  should  be  consistent  with  the  system 
being  studied.  If  U(r)  is  to  be  a liquid  metal  pair  potential  then 
F (q)  should  be  calculated  for  the  appropriate  liquid  metal  density. 
Pair  potentials  for  Li,  Na,  and  K were  calculated  using  Equation  (34). 
The  results  are  given  in  Figure  11.  The  F(q)  used  to  obtain  these 


U(r)  were  calculated  with  the  same  treatment  of  the  conduction-core 


vr» 
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exchange  and  dielectric  screening  that  worked  successfully  in  the 
solid  metal.  However,  for  the  liquid  calculation  F(q)  used  a volume 
per  ion  that  was  the  same  as  the  observed  experimental  density  of  the 
liquid . 

The  U(r)  curves  in  Figure  11  show  several  qualitative  features 
that  should  be  noted.  First,  there  is  a strongly  repulsive  core  region 
that  will  keep  pairs  of  ions  from  coming  too  close  together.  Also, 
there  is  a strong  potential  minimum  which  occurs  at  a radius  that  is 
in  the  vicinity  of  the  nearest  neighbor  or  first  coordination  shell 
distance.  If  the  liquid  were  arranged  in  a crystalline  structure, 
the  minimum  of  the  U(r)  in  Figure  11  occurs  at  a distance  that  lies 
between  the  first  and  second  nearest  neighbor  distances.  Finally,  at 
large  r,  U(r)  exhibits  an  oscillatory  behavior  which  is  a characteristic 
of  metallic  pair  potentials  and  is  due  to  the  screening  properties  of 
the  conduction  electron  gas  (March,  1968) . 


* 
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IX.  MONTE  CARLO  RESULTS  FOR  THE  LIQUID  STRUCTURE 

The  U(r)'s  in  Figure  11  were  used  to  determine  the  cummulative 
distribution  function  G(r)[i.e.  equation  (33)1  using  a Monte  Carlo  calcu- 
lation. These  functions  for  liquid  Li,  Na  and  K were  computed  at  the 
melting  temperatures  and  corresponding  densities  (see  Table  3) . The 
Monte  Carlo  calculations  were  done  for  systems  of  125  and  216  particles 
initially  in  a cubic  configuration.  The  first  10,000  configurations 
were  run  without  averaging  to  let  the  particles  evolve  to  more  real- 
istic configurations.  Then,  an  additional  50,000  to  100,000  configura- 
tions were  run  to  compute  the  ensemble  average  G(r)  functions  for  Li, 

Na  and  K.  Due  to  the  use  of  the  minimum  image  convention,  G (r)  is 
defined  over  the  finite  range  0 -*■  L/2  (Wood  and  Parker,  1957). 

Equation  (35) 

g(r)  = l/(4TTr2n0)8G(r)/3r  (35) 

was  used  to  calculate  g(r)  from  G(r)  and  the  results  of  these  calcula- 
tions for  Li,  Na  and  K are  displayed  in  Figures  12,  13  and  14  respec- 
tively. These  curves  represent  the  216  particle  system.  Results  for 
the  125  particle  calculations  were  not  significantly  different/  indicat- 
ing that  many-body  effects  are  adequately  accounted  for.  It  is  to  be 
noted  that  the  radial  distribution  curves  exhibit  characteristic 
behavior,  that  is,  at  small  r,  g(r)  is  zero  due  to  core  repulsion.  A 
sharp  peak  occurs  in  g(r)  at  intermediate  r and  corresponds  to  a 
coordination  shell  of  nearest  neighbors.  At  large  r values  (i.e., 
beyond  the  first  peak)/  g(r)  approaches  the  ideal  liquid  limit  of  1 
since  at  these  distances  the  liquid  begins  to  look  more  and  more  like 


a continuous  isotropic  medium. 


function  for  Li.  The  curve  is  the  result  of  a Monte 
ing  the  pair  potential  in  Figure  11. 
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Table  3.  Temperature  and  density 
of  Li,  Na,  and  K at 
their  melting  point. 


T ( °K) 
m 

nQ (gm/cm3 ) 

Li 

4 54 

0.515 

Na 

373 

0.928 

K 

337 

0.828 

The  g(r)  curves  have  been  extended  beyond  the  limits  of  our  Monte  Carlo 
results  (solid  part  of  the  curve)  through  the  use  of  the  Percus-Yevick 
theory  which  predicts  that  at  large  r,  g(r)  should  behave  like,  March 
(1968) 

g ( r)  = 1 + Bcos  (2Kf r + 4>)/r3  , (36) 

where  B and  ip  were  chosen  so  that  there  was  a smooth  transition  between 
the  Monte  Carlo  and  analytic  protions  (dotted  part)  of  the  g(r)  curves. 

The  structure  factor,  S(q),  was  calculated  for  the  simple  alkalis 
using  Equations  (26)  and  (27).  The  S(q)  curves  in  Figures  15,  16,  and  17 
for  Li,  Na  and  K, respectively,  were  obtained  for  q > 2K^  by  Fourier 
transforming  g(r)-l.  At  low  values  of  q [i.e.,  below  the  first  peak 
in  S(q)],the  Fourier  transform  method  is  not  accurate.  Following  the 
suggestion  of  Fowler  (1973), S(q)  was  calculated  at  low  q directly  from  its 
definition  in  Equation  (26) . The  explicit  expression  used  to  evaluate 
S (q)  at  low  q was 

3 N 

S (q)  = 1/3  l l (sin2q  x + cos ‘q  x .) 


(37) 
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where  a = 1,  2,  3 denote  x,  y,  and  z components,  i is  a particle  label, 
and  q was  chosen  as  an  integer  multiple  of  2ti  (i  + j + £)/L.  Note  that 
this  is  along  a cube  diagonal.  As  before,  the  bar  denotes  an  ensemble 
average  which  was  computed  by  using  the  Monte  Carlo  technique.  The  use 
of  Equation  (37)  at  low  q to  calculate  S (q)  is  much  more  accurate  than 
the  Fourier  transform  method.  However,  if  S(q)  is  calculated  from 
Equation  (37) , there  is  a problem  with  resolution  since  the  q values 
must  be  interger  multiples  of  2TT/L  due  to  the  periodic  boundary  conditions. 
Calculating  S(q)  at  low  q from  Equation  (37)  is  more  accurate  in  the 
region  0 < q < 2K^.  , but  this  method  will  provide  only  a limited  number 
of  points  with  which  we  can  define  the  S(q)  curve.  Table  4 gives  the 
Monte  Carlo  low  q results  for  the  S (q)  of  Li,  Na , and  K.  The  experimental 
points  in  the  table  are  interpolated  from  the  experiments  of  Greenfield 
et  al.  (1971) . 

The  S(q)  values  for  Na  and  K are  in  reasonable  agreement  with 
experimental  values  as  can  be  seen  from  Figures  16,  16,  and  17.  At  low 
q,  where  S(q)  becomes  quite  small,  the  percentage  error  is  relatively 
larger  but  the  qualitative  agreement  with  experimental  data  is  still 
excellent.  S(q)  should  approach  the  compressibility  limit  at  low  q 

S (q)  * noK0TKt  as  q ->  0 , (38) 

where  is  the  isothermal  compressibility,  T is  the  temperature,  and 
K0  is  the  Boltzmann  constant.  Greenfield  et  al . list  the  q -*  0 limits 
of  S (q)  for  Na  and  K as  0.0240  and  0.0247,  respectively.  Our  calculated 
S(q)'s,  though  not  in  precise  agreement  with  these  limits,  are  still  more 
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Table  4.  S(q)  values  for  bi,  Nu,  and  K at  low  q as  predicted  by  Monte 
Carlo  computations.  Units  of  q are  inverse  Bohr  units. 


Li 

Na 

K 

q 

S (q) 

q 

S fq il 

Exp. 

q 

S(q) 

Exp. 

0.  196 

0.0225 

0.160 

0.0318 

0.0263 

0.129 

0.0212 

0.0286 

0.393 

0.0262 

0.321 

0.0324 

0.0321 

0.258 

0.0216 

0.0353 

0.  589 

0.0291 

0.481 

0.0415 

0.0452 

0.387 

0.0266 

0.0479 

0.  786 

0.0703 

0.642 

0.0750 

0.0800 

0.517 

0.0492 

0.0800 

0.983 

0.246 

0 . 802 

0.152 

0.  196 

0.646 

0.133 

0.218 

0 . 96  3 

1.07 

1.07 

0.775 

1.17 

1.25 

accurate  than  results  employ inq  the  transform  method  to  calculate  S(q) 
(Murphy,  1977;  Murphy  and  Klein  L973) . The  analytic  extension  of  the 
q(r)  curves  heiqhtens  the  first  peak  of  the  S(q)  curve  and  lowers  the 
first  minimum  slightly  but  does  not  significantly  affect  the  rest  of 
the  S(q)  curve.  Hence  the  extension  of  g(r)  makes  some  quanitative 


improvements  in  S(q)  but  does  not  make  any  qualitative  changes. 
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X.  LIQUID  METAL  ELECTRICAL  RESISTIVITY 


The  electrical  resistivity  p for  Li,  Na  and  K has  been  calculated 

Li 

using  Ziman's  theory  (Ziman,  1972) 

2kc 


i - 3"  i" 

l 4tle‘vf‘n0k{. ' 


S(q)  |«(q)  I ‘ q 3 dq 


(39) 


The  equation  is  derived  in  Appendix  E where  the  symbols  are  defined.  The 

input  required  for  this  calculation  is  the  structure  factor  S (q)  and  the 

-+  ♦ . . -*• 

pseudopotential  matrix  element  <k  + q|w|k>  evaluated  on  the  Fermi  sphere 

(i.e.,|k|  = |k  + q|  = kf)  . Table  5 gives  the  calculated  and  experimental 

values  of  p for  the  simple  alkalis.  The  computed  values  of  P.  were 
L Lj 

obtained  with  our  theoretical ly  determined  value  of  S(q).  For  the  fol- 
lowing reason,  this  constitutes  a more  rigorous  test  of  the  pseudopotential 
than  if  the  experimental  S(q)  was  used.  S(q)  was  calculated  using  the  same 
pseudopotential  as  that  used  to  describe  the  scattering  potential  in  the 
calculation  of  the  resistivity.  This  procedure  gives  us  a degree  of 

"self-consistency"  vis-a-vis  the  ion-ion  potential  in  the  resistivity 

th  exp 

calculations.  The  differences  between  p and  p are  not  unreasonable 

L Li 

when  one  considers  the  extreme  sensitivity  of  P^  to  both  S(q)  and  the 
matrix  element  |'.kjw|k  > q>.  Our  calculations  of  p for  various  nvitrix 
elements  have  stiown  up  to  factor  of  two  change  in  P , though  the  matrix 

Lt 

elements  at  q = 0 were  of  the  order  of  -0.1  rydbergs  and  -0.01  rydbergs 

at  q = q . The  differences  were  less  than  1*  at  low  q and  never  exceeded 
kf 

0.01  rydbergs  in  the  2k^.  interval  for  q,  yet  the  computed  resistivities 
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Table  5. 

Theoretical  and  experimental 
values  for  the  electrical 
resistivity  of  Li,  Na,  and  K 
at  their  melting  points. 

exp 

hL 

(Uft-cm) 

(pft-cm) 

Li 

24.0 

16.  5 

Na 

9.7 

o 

CO 

K 

13.0 

8.2 

differed  by  a factor  of  two.  The  reason  for  this  extreme  sensitivity  is 
that  the  resistivity  depends  on  the  product  of  S (q)  times  |<k  + q|w|k>| 
in  the  region  0 i q < > kf.  S(q)  is  sharply  peaked  right  around  q = 2kf 
and  will  act  something  like  a delta  function.  However)  the  matrix  element 
undergoes  a sign  change  in  the  2k^  region  and  is  rather  small.  Con- 
sequently, if  two  matrix  elements  differ  by  a small  absolute  amount  in 
this  region,  the  resistivities  can  be  significantly  different. 

The  matrix  element  in  Equation  (39)  assumes  a rather  simple  form 
because  it  is  evaluated  on  the  Fermi  sphere  [<k  + q | w|  K>-*-w  (q)  ] . It  is 
tempting  to  evaluate  W (q)  using  the  local  approximation  (ignore  k dependence 
of  the  matrix  element)  since  we  are  constrained  to  the  Fermi  surface.  How- 
ever, this  local  approximation  significantly  affects  the  evaluation  of  the 

s c 

screening  potential  W (q) , which  is  itself  local  but  depends  on  the 
nonlocal  parts  of  the  pseudopotential  [see  Equation  (D.ll)].  If  the  k 
dependence  of  the  bare  matrix  element  is  ignored  in  evaluating  the 
screening  field,  it  can  produce  up  to  40%  differences  in  the  calculation 
of  p . Even  for  a simple  metal  like  Na  we  find  that  the  local  approxima- 

li 

tion  leads  to  significant  errors  ( > 20%)  in  calculations  of  the  phonon 
spectra  for  the  solid  metal.  Therefore  the  nonlocality  of  the  pseudo- 


potential  was  retained  in  our  calculations  of  the  resistivity. 
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XI.  POSSIBLE  IMPROVEMENTS 

Consistently  good  results  have  been  obtained  for  elastic  shear 
constants,  phonon  spectra,  and  liquid  metal  structure  of  the  simple 
alkali  metals  using  a fully  non  local  first  principles  pseudopotential 
formalism.  During  the  course  of  this  investigation  procedures  for 
improving  the  pseudopotential  calculations  have  suggested  themselves. 

In  particular,  a possible  means  of  improving  the  core  electron  states 
is  considered. 

A first  principles  pseudopotential  is  determined  to  a large  extent 
by  the  choice  of  wave  functions  for  an  ion  in  free  space.  This  is  a 
reasonable  approximation  for  the  alkalis  but  one  must  wonder  if  a better 
approximation  could  not  be  found.  Sometimes  the  free  atom  wave  functions 
are  used  but  these  do  not,  in  principle,  constitute  any  fundamental 
improvement  as  a representation  of  the  core  electron  states  in  the 
crystal  and  gave  greater  difficulty  when  calculating  core  shifts.  We 
would  now  like  to  suggest  as  an  improved  representation  of  the  core 
states  in  a metal  the  so-called  pseudo-atom. ' A 'pseudo-atom'  is  an 
ion  which  is  surrounded  by  a compensating  screening  charge  cloud  such 
that  the  ion  and  screening  charge  constitute  a neutral  charge  system. 

This  is  approximately  the  case  in  a metal  where  the  conduction  electrons 
screen  the  ionic  potentials.  The  well-known  Herman-Ski liman  (1963)  atomic 
structure  program  can  be  easily  modified  so  that  pseudo  atom  wave 


functions  can  be  determined  self  consistently. 
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J.  A.  Moriarty  (1974)  has  calculated  and  used  simple  pseudo-atom 
wave  functions  successfully  in  pseudopotential  calculations.  Moriarty 's 
pseudo-atom  was  more  simplistic  than  what  we  are  proposing.  In  his 
calculations*  the  screening  charge  of  conduction  electrons  was  approximated 
by  a rigid  sphere  of  uniform  charge  density.  The  wave  functions  for  an  ion 
centered  in  this  sphere  were  then  calculated  self-consistently  using  a 
modified  Herman-Skil lman  program.  One  can  go  a step  further,  however, 
and  calculate  the  wave  functions  of  an  ion  surrounded  by  a spherically 

SC 

symmetric  screening  charge,  where  the  screening  potential  V can  be 
expressed  in  terms  of  the  bare  ion  potential  and  an  appropriate  dielectric 
function.  It  will  presently  be  shown  how  this  may  be  accomplished. 

First  note,  however, that  the  self-consistently  calculated  pseudo-atom 
wave  functions  can  be  applied  to  more  than  just  first  principles 
pseudopotential  calculations.  Improved  core  electron  states  can  be 
used  in  APW  and  tight  binding  band  structure  calculations  Calloway 
(1976).  Also, a self  consistent  pseudo-atom  should  make  an  excellent 
representation  for  an  impurity  in  a metal. 

The  Herman-Skil lman  program  is  designed  to  solve  the  radial  part  of 

the  Schrodinger  equation  self  consistently  for  the  wave  function  P „(r) 

nx. 

P „ (r)  = E P (r)  . (40) 

nx.  nx  n£ 

Equation  (40)  is  the  radial  Schrodinger  equation,  is  the  eigenvalue 

for  P . (r),  and  the  potential  v(r)  is  normally 
nX 


, tlSUL  ♦ v.(r, 

r2 


r 
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vtr)  - ^ ♦ 2f  "<^d,r  ♦ vex(r) 
r J I r-r ' | 


(41) 


The  first  term  in  Equation  (41)  is  the  nuclear  coulomb  potential,  and 
the  second  term  is  the  coulomb  potential  of  the  spherically  averaged 
core  electron  charge  density 


occup. 

nlt’  ■ V p^<r)/4’t2 

nx, 


(42) 


where  M ^ is  the  occupation  number  for  orbital  n,  &.  The  third  term  is 

the  exchange  potential  of  the  electrons  which  is  normally  approximated 

by  a local  potential.  Equation  (40)  can  be  solved  in  an  iterative 

manner  for  the  P „(r^  and  E .of  a free  ion  or  atom.  The  Herman- 
n,Jt  n,x, 

Skillman  program  can  be  converted  to  a pseudo- atom  calculation  by  replac- 
ing V°(r)  by  V(r)  in  Equation  (40).  The  pseudo-atom  potential  V(r)  is 


V(r)  = vu(r)  + vsc(r)  , 


(43) 


composed  of  the  bare  potential  V°(r)  from  Equation  (41)  plus  a screening 

sc  sc 

potential  V (r) . V (r)  may  be  obtained  from  V(r)  in  the  following 
manner.  First,  Fourier  transform  V(r)  to  obtain  V(q) . 


V u(q)  = 


V“(r)  e '*''1  rdJr  . 


(44) 


SC  o 

V (q)  can  be  expressed  in  terms  of  V (q)  as 


V (q)  = V°  (q) [ 1-C* ( q ) ] /e* (q)  , 


(45) 


- A 
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where  e (q)  [see  Equation  (D.12)]  is  a static  dielectric  response  function 

which,  in  our  calculations,  was  chosen  as  the  SSTL  dielectric  function. 

s c 

V (r)  in  v quation  (43)  is  obtained  from 


sc 

V (r)  = 


(271) 


VSC(q)elq,rd3q 


(46) 


which  is  simply  the  inverse  transform  of  Equation  (44) . The  modifications 
to  the  Herman-Skillman  program  allow  one  to  calculate  the  self  consist- 
ently screened  wave  functions  of  the  pseudo-atom  in  the  linear  approxi- 
mation. 


The  pseudo-atom  wave  functions  for  Li  were  calculated  to  test  the 
feasibility  of  this  method.  As  one  would  expect,  the  first  eigenvalue 
for  the  pseudo-atom  was  significantly  less  in  magnitude  than  the  ionic 
eigenvalue.  This  is,  of  course,  due  to  the  screening  potential  in  the 
pseudo-atom  calculation  which  causes  the  core  charge  to  relax.  It  would 
be  interesting  to  make  a systematic  study  of  the  use  of  pseudo-atom 
wave  functions  in  first  principles  pseudopotential  calculations.  (There 
are  advantages  such  as  the  following.)  A 'core  shift'  is  calculated 
for  the  ionic  eigenvalue  that  is  the  interactions  between  the  core 
charge  and  the  10PW  conduction  charge  density.  In  both  Moriarty.4s  and 
our  pseudo-atom  calculations,  this  interaction  is  taken  into  account  in 
the  wave  function  calculation,  in  a self  consistent  manner;  thus,  the  use 
of  the  pseudo-atom  leads  to  the  elimination  of  the  core  shift  in  pseudo- 
potential calculations.  Our  pseudo-atom  calculation  would  not  be 
consistent  with  the  sphere  model,  but  it  is  hoped  that  the  inclusion  of 
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screening  in  the  wave  function  calculation  is  a more  realistic  approach. 
Finally, we  note  that  one  could  in  principle  achieve  a degree  of  self 
consistency  between  the  core  wave  functions  and  the  crystal  potential 
by  performing  an  iterative  type  of  calculation,  where , at  each  iteration, 
the  core  states  are  calculated  in  the  field  of  crystalline  potential 
(spherically  symmetric  part  only)  from  the  previous  iteration.  This  would 
probably  require  a very  large  amount  of  computation  but  could  be  done. 
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XII.  SUMMARY  AND  CONCLUSION 

Using  a first  principles  fully  nonlocal  pseudopotential,  a systematic 
study  of  the  simple  alkali  metals  has  been  made.  The  phonon  spectra 
predicted  by  the  pseudopotential  were  shown  to  be  sensitive  to  the  exact 
forms  of  the  core  electron  states  and  crystalline  potential  used  to 
construct  the  pseudopotential.  In  particular,  it  was  demonstrated  that 
standard  approximations  used  to  describe  the  orthogonalization  hole 
contribution  to  the  conduction  electron  potential  are,  in  general,  in- 
accurate. The  exact  treatment  of  the  orthogonalization  hole  potential 
was  shown  to  produce  significant  changes  in  calculated  lattice  dynamical 
properties  of  the  solid  alkali  metals  Li,  Na,  and  K. 

Computer  simulations  of  the  alkalis  in  the  liquid  state  were  per- 
formed using  pair  potentials  derived  from  our  HFP  fully  nonlocal 
pseudopotential  calculations.  Radial  distribution  functions,  g(r),  and 
the  static  structure  factor,  S (q) , were  computed,  using  a Monte  Carlo 
technique.  Very  good  agreement  between  experimental  and  theoretical 
S (q)  curves  was  obtained.  Liquid  resistivity  calculations  were  per- 
formed with  a consistent  treatment  of  the  structure  factor  and  scat- 
tering matrix  elements.  Specifically,  Equation  (39)  utilized  a 
structure  factor  S (q)  that  was  predicted  by  the  same  pseudopotential 
matrix  element  W(q)  used  to  describe  the  electron  scattering  in  the 
liquid.  Since  the  completion  of  this  thesis  research,  Sun  (1978) 
has  shown  that  additional  refinements  to  the  present  work  can  further 
improve  the  agreement  between  the  calculated  transport  properties, 
resistivity,  thermoelectric  power,  thermal  conductivity,  transport 


parameters,  and  experiment. 
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Finally,  suggestions  were  made  for  improvements  in  calculations 
of  the  core  electron  states.  A method  for  calculating  self  consistently 
screened  pseudo-atom  wave  functions  was  derived.  These  wave  functions 
should  constitute  a better  representation  of  the  core  electron  states 
in  a crystal. 
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APPENDIX  A 

THEORY  OF  PHONCNS  FOR  RCC  STRUCTURE 

In  this  appendix*  it  will  be  demonstrated  how  the  phonon  frequen- 
cies depend  upon  the  enemy  wave  number  characteristic  F(|q|).  The 
fol  lowing  discussion  is  limited  to  the  case  of  materials  with  one 
atom  per  unit  cell  but  may  easily  be  qeneralized  to  any  number  of 
atoms  per  unit  cell.  In  the  harmonic  approximation  the  potential 
energy  of  a crystal  with  one  atom  per  unit  cell  may  be  written  as, 


* = + ~ y )'  'J'  n (O  u («  ')  , 

0 '■  ‘i  , (Xp  (X  Is 


?, , a < B 


(A.  1) 


where  '1  is  explicitly 

an 


VM'> 


,vv 

aua(f)  9U6(JI') 


(A.  2) 


a(5) =0 
Up(t')»0 

The  subscripts  oi  and  B denote  spatial  directions,  u^^  and  llp(C'l 

are  the  displacements  of  atoms  C and  C ' from  equilibrium,  and  the 

c 

second  derivative  of  the  crystalline  potential  V is  evaluated  at 

-> 

equilibrium.  An  equation  of  motion  for  the  displacements  u ^ can  be 
obtained  by  taking  the  qradient  of  Equation  (A.l) 

" I * a,.  (A-3) 


mug(«)  ’ " /•  afl<«,0  UB(«') 

V t P 


• » >n  sp«»et  ra  are  obtained  by  calculat  ing  t he  normal  mode  solu- 

-*■ 

. Equation  (A.l);  for  a phonon  of  wavevector  Q and  frequency 
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1 „ ,*x  -i«t 

uix((i)  ‘ / Ua(Q)  e 

HU 


where  x^  is  a lattice  vector.  Equation  (A.  3)  can  be  interpreted  as 
representing  the  force  on  atom  i in  the  a'th  direction  due  to  dis- 
placing the  atoms  i ' in  the  j3  direction.  Substituting  Equation  (A. 4) 
into  Equation  (A. 3)  yields  the  well  known  secular  equation 


W2(y)  Uo(Q)  - l DuB(Q)  Ug(B)  . 


The  dynamical  matrix  elements  D (Q)  are  defined  as 


“os*  W1-*’’ 


-iQ* (x£-xe,) 


The  phonon  frequencies  t o(Q)  are  computed  by  evaluating  the  secular 
determinant  that  is  implicitly  contained  in  Equation  (A. 5).  The  rest 
of  this  appendix  will  show  how  the  dynamical  matrix  is  evaluated. 

The  elements  of  the  dynamical  matrix  have  contributions  from  two 

different  sources.  These  are  the  electrostatic  energy  of  the  ions  in 

-► 

the  crystal  and  the  energy  of  the  conduction-electron  gas.  P is 

CXP 

o s ce  ► 

the  sum  of  these  contributions  which  are  labelled  P , (0)  and  D ,,  (O)  , 

txp  exp 

respectively.  The  electrostatic  contribution  will  be  discussed  first. 


ELECTROS TAT 1 C CONTRIBUTION  TO  P AQ) 

txp 


The  charge  centered  on  the  lattice  site  x „ , will  interact  with 


the  point-ion  in  the  i'th  unit  cell.  To  compute  the  energy  of  this 
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interaction* we  must  know  the  Coulomb  potential,  v (x) , associated  with 

->  gg  -> 

the  charges  centered  at  . The  potential  v (x)  is  due  to  the 

charge  density 

p(x)  = z *e  [6  (x-x^,  , ) - 1A2]  , (A.  7) 

which  is  composed  of  a positive  point  ion  of  magnitude  z*e  and  a 

compensating  uniform  electron  charge  density.  A gaussian  charge 

density  of  half  width  n is  added  and  subtracted  from  Equation  (A. 7)  to 

make  lattice  sums  convergent.  This  procedure  does  not  alter  any  of 

► 

the  physical  results.  The  charge  density  p(x)  then  becomes 


p(x)  = z e 


3 3/ 2 

n IT 


exp[  (-x-Xjj , )  l  2/n2 


1/SI 


<5  (x-xD , ) 


,3n  V2 


exp  l(x-xr)2/n 


= Pj  (x)  + p2  (x)  , (A. 8) 

-►  -> 

where  pj (x)  and  P2 (x)  are  the  charge  densities  in  the  first  and  second 

l 1 brackets,  respectively . The  Coulomb  potentials,  v ( (x)  and  v2  (x) , 
->  > 

associated  with  pj (x)  and  p. (x)  are  evaluated  as  follows.  The  distri- 
bution p, (x)  is  spherically  symmetric  and  Gauss'  law  may  be  used  to 

find  its  electric  field.  Integrating  over  the  electric  field  yields 
the  potential 


es , »•. 
v,  (x) 


* 

2 e 


erfcjx-x, , | /n 


x-x„ 


(A. 9) 
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Equation  (A. 9)  excludes  a constant  term  which  disappears  upon  differen- 
tiation (Tosi,  1964). 


It  is  convenient  to  express  pt (x) 


Pt  (x) 


2*e  -iq*(xr-x) 

L ~cT  e 
q^O  * 


as  a Fourier  series 

2 2 -*  -*• 

-q  n /4  iq*x 
e e 


(A. 10) 


where  the  q = 0 term  in  Equation  (A. 10)  is  missing  because  p( (x' 
represents  a neutral  charge  system.  Using  Poisson's  equation 


es  , 

v (q) 


- 4ir  p.  (q)/q‘ 


es , 


one  obtains  for  v (x) , 


(A. 11) 


es 

v (x) 


l 

q*0 


4iiz*e 

q‘il 


-qV/4 


-iq* (x^,-x) 


* 

z e 


;rfc 


V 


/n 


V 


(A. 12) 


The  dynamical  matrix  elements  in  Equation  (A. 6)  may  be  rewritten 


as 


GS 

.es  1 v d V (x) 

aP  Q = m ^x^3xp 


4 4 4 

■iQ* (x-x^,) 


(A. 13) 


X=x, 


GS  4 

In  Equation  (A. 13) , the  derivatives  of  V ‘ (x)  are  taken  with  respect 

4 

to  x and  not  the  displacements  u as  in  Equation  (A. 6).  The  function 

gs  4 es  4-  . 

V (x)  differs  from  v (x)  by  a factor  z e.  Hence,  by  substituting 

z*e  times  Equation  (A. 12)  into  Equation  (A. 13)  one  obtains 
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»3®  ■si.'iJ 


4Tlz*2e2 


ft'  q^O  q m 


-qV/4  •iq,(x£'“x£) 

q q ,,  e e 

a 3 


* 2 2 

- z e 


^2  erfc|x-Xj, |/n 


9x  9x0 

a g 


Ix-X^, 


> e 


-iQ* <x^"x£') 


-»  > 
x=*x„ 


(A. 14) 


It  is  apparent  that  the  above  expression  for  is  a function  of 

x^-x^,  which  must  also  be  a lattice  vector.  Therefore,  in  Equation 

**►  r»  r* 

(A. 14),  it  is  possible  to  replace  x.-x^,  by  x.  and  l by  l . Note  that 
when  £ = £',  the  second  term  in  Equation  (A. 14)  is  singular.  This 
difficulty  may  be  circumvented  by  making  use  of  the  condition 


V*'r) 


l 

I'/l 


W4'1') 


(A. 15) 


which  results  from  the  fact  that  the  system  is  invariant  under  a 
uniform  translation.  Utilizing  Equation  (A. 15)  and  the  above  con- 


siderations results  in  Equation  (A. 16)  for  D 


es 

oi3 


es 

ci3 


.2„2 


is-  i j q e q n /4  lq  **•  + _il ££.fc_  ixI/n 

” t-o  i Jo  q'n % 11  Vxs  III 


1 - e 


-iQ*x£ 


(A. 16) 


Equation  (A. 16)  may  be  simplified  by  using  the  lattice  sum  rule 


N -iq*x. 

} e ■ N 6 

£=0 


q,G 


(A. 17) 


r 
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Making  use  of  Equation  (A. 17)  to  eliminate  the  £ sum  in  the  first  term 
of  Equation  (A. 16),  and  explicitly  differentiating  the  erfc,  one  obtains 


es  -► 

DaB(Q)  = mrt 


47TZ*2e2  y' 


(G^)C,<G^’b  „-|&3|V/4  _ Vb  e-G2n2/4 
|g+q|  g2 


z*v  y < 

I x£,ax£, B 

3.rfo|;,|/n  -*l/n  -Xt/r|Z 

1 £ 6e  4e 

m Vo  | 

, ’i 

+ + 

x^  n/n  x2  n3/ir 

- 6 


aB 


2e  Xjl/n  erfc|x£|/n 


nx 


.3 


iQ*x0 


1 - e 


(A. 18) 


The  prime  on  the  G sum  indicates  that,  when  G = 0,  the  term  with  1/G2 
is  eliminated. 


THE  CONDUCTION  ELECTRON  CONTRIBUTION  TO  D 


aB 


The  conduction  electron  contribution  to  D „ will  now  be  considered. 

aB 


The  energy  of  the  conduction  electron  gas  depends  on  the  crystal  struc- 
ture through  Equation  (18)  which  we  duplicate  here  for  convenience 


E = l S (q)  S*  (q)  F (q) 

q 


(A. 19) 


The  structure  factor  S,  . is  a function  of  the  lattice  coordinates. 

(q) 


For  a crystal  with  one  atom  per  unit  cell,  it  may  be  written  as 


- lr 

s q =n  \e 


(A. 20) 


1 


1 
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where  N is  the  number  of  unit  cells,  x^  is  the  position  of  atom  i,  and 
u^  is  its  displacement  from  equilibrium.  It  is  convenient  to  expand 
S(q)  in  terms  of  the  displacement,  u^ 


1 r iq’xH 

s (q>  = n l e 

N i 


i - icru£  + 


(-iq*u^) 

21 


(A. 21) 


For  a pure  normal  mode  of  wave  vector  Q,  the  lattice  displacement  u^ 
may  be  written  as 


u£  = 


U(Q)  e + U (Q)  e 


(A. 22) 


In  Equation  (A. 22)  the  normal  mode  displacement  has  been  expressed  in 
a form  that  ensures  that  it  is  real.  Substituting  u^  into  Equation 
(A. 21)  and  using  Equation  (A. 27),  S( q)  becomes 


S (q) 


U<2)  <5+  + - 

q-Q,G 


-iu)t 

e 


-*•*  •> 

+ u (Q) 


q+Q,G 


iu)t 

e 


-1 

2m 


[q*U (Q) ] 2 6 


-*■  -+• 


-2ia)t 
e + 


2q*U 


q*U 


+ • • * . (A. 23) 

The  product  S* (q)  in  Equation  (A. 19)  will  be  nonzero  only  when  the 
Kronecker  deltas  in  the  product  have  similar  arguments.  Retaining 
only  terms  up  to  second  order  in  the  displacements,  one  obtains  the 
following  for  the  conduction  electron  energy  change  due  to  the  normal 
mode  displacements 


+ tq 


■u  (Q)]2 


,2icot 


q-2Q,G 


1 


ae  = - 


“ JT  I l (G+g)  (G+Q)g  U*(Q)  U (Q)  F(|g+q|) 

3 la, 6 a B a p 

+ I (G-Q)0(G-C)g  Ua(C)  U*(C)  F(|g-q|) 

a,  3 


" 2 I Ga  Gg  Ua(Q)  UJ(Q)  F(|G|) 


(A. 24) 


In  Equation  (A. 24) , a and  3 denote  the  cartesian  components  x,  y,  z, 

“►  ■+■  •>  -f 

and  the  factors  G + Q and  G - Q are  equivalent  because  of  the  sum  over 

G-vectors.  Equation  (A. 24)  represents  the  energy  difference  of  the 

, -► 

conduction  electrons  between  the  equilibrium  lattice  (u^  = 0)  and  the 

distorted  lattice.  An  equivalent  expression  in  terms  of  the  dynamical 
matrix  can  be  obtained  by  combining  Equations  (A.l) , (A. 6) , and  (A. 22) 
(include  a sum  over  Q)  to  find  the  second  order  energy  change  per  ion 
caused  by  a lattice  displacement.  This  is  given  by 


* - <t>°  k l Dr»R  (®)  + D*g(Q>  Ua(Q)  U*(Q)  . (A.  25) 


2 a3  aS 


Comparing  Equations  (A. 24)  and  (A. 25)  results  in  the  conduction  electron 
contribution  to  the  dynamical  matrix.  This  is  given  by 


Da3(Q)  = m E (G+Q)„(G+Q)ft  F(|g+q|)  - G^Gft  F(|g|)  , (A. 26) 


'a' 


->  I -►  I 

where  the  G = 0 term  involving  F(|G|)  is  omitted  from  the  sum.  Hence, 


with  Equations  (A. 26)  and  (A. 18) , the  phonon  frequencies  u>{Q)  can  be 
determined  by  evaluating  the  secular  determinant  in  Equation  (A. 5) . 
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APPENDIX  B 


ELASTIC  CONSTANTS  FOR  BODY  CENTERED  CUBIC  MATERIALS 


In  general  the  elastic  energy  density  W is  related  to  the  elastic 
constants  CL  ^ and  the  strains  e^  through  the  following  relationship 


w “ W.  + =-  y C . . e . e . 
° 2 i,j=l  13  1 3 


This  equation  is  only  valid  for  small  deformations  and  is  bar  3 on 
Hooke's  Law  or,  equivalently,  the  harmonic  approximation.  W0  repre- 
sents the  energy  per  ion  of  the  undeformed  (i.e.,  equilibrium)  crystal. 
The  strains  e.  are  defined  as  follows 

l 


3x  * 


3v  9w  _ 3w  9v 

3y  ' 6 3 3z  ' 64  3y  + 3z  ' 


3u  3w  3u  3v 

es  " 3z  + 3x  and  es  = 3y  3x  ' 


where  u,  v,  w are  the  x,  y,  and  z components  of  the  displacement  in  a 
crystal  caused  by  a deformation,  , e2 , e3  are  related  to  distortions 
of  the  x,  y,  and  z axis,  and  e^ , e5,  e^  describe  angular  deformations  of 
the  angle  between  the  axis  defining  the  yz , xz , and  xy  planes,  respectively. 

It  can  be  shown  that,  at  most,  there  are  21  independent  elastic 
constants.  Body  centered  cubic  (bcc)  materials  have  only  three  independent 
elastic  constants.  Thus, for  a bcc  crystal,  it  suffices  to  calculate 
three  constants  designated  by  Cn,  C12,  and  C41(  in  order  to  specify 
the  elastic  properties  of  the  material  since  the  other  elastic  con- 
stants are  either  zero  or  equal  to  one  of  these.  In  the  rest  of  this 
section,  two  specific  volume  conserving  strains  will  be  considered 
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that  can  be  related  to  CM~C12'  and  c,,4  with  the  aid  of  Equation 
(B.l)  . 

Consider  the  vectors  describing  a primitive  unit  cell  for  a bcc 
crystal  before  and  after  a deformation.  If  a1#  a2,  a3  are  the 
primitive  vectors  before  the  strain,  and  a , a2,  a3  are  the  vectors 
after  the  strain,  then  they  are  related  by  the  strain  matrix  e 


" ai 

s 

= e 

ii 

t «T 

(UI' 

Vl' 

Wl) 

-4- 

" a2 

= e 

-► 

a2  = 

(u2 , 

V2' 

W2> 

■*/ 
a 3 

-4 

- 3 3 

= e 

3 3 “ 

(U3' 

V3' 

w3) 

(B.3) 


where  the  strain  matrix  e is 


e = 


®1 

e6/2 

e5/f 

e6/2 

e2 

e4/2 

• 

e5/2 

e4/2 

e3  _ 

(B.4) 

Hence,  if  the  vectors  a.  and  af  are  known,  then  the  strains  e,  can  be 

' 1 i ' l 

obtained  through  Equations  (3)  and  (4)  as  there  are  nine  equations  and 
only  six  unknowns. 

The  bcc  primitive  vectors  in  the  unstrained  state  are 

ax  - | (1,1,-1)  , a2  = | (-1,1,1)  and  a3  = | (1,-1,1)  » (B.5) 

where  a is  the  lattice  constant.  Now,  if  the  z axis  is  stretched  and 


the  x and  y axis  compressed  so  that  the  primitive  cell  volume  is 


~r-  ~r-  -T  I I f 

conserved  (see  Figure  18)  (i.e.,  a1*a2  x a3  = a^a,  * a3),  then  the 

vectors  may  be  expressed  as 


a . i/2  i/2  -i. 

a 1 = 2 (ei  '£i  '-£i  > 


a , i/2  1/2  -d 

J ("ei  «ei  'ei  > 


a . 1/2  i/2  -d, 

7 '~e,  -e.  > 


Ej  is  the  strain  parameter  which  is  equal  to  1 at  equilibrium  and  i is 
arbitrary  and  can  be  chosen  for  convenience.  Substituting  Equations 
(B.4) , (B.5),  and  (B.6)  into  Equation  (B.3)  results  in  a system  of 
coupled  equations  for  the  e^.  Solving  this  set  of  simultaneous  equa- 
tions, one  obtains  the  strain  matrix  elements  in  terms  of  the 
parameter  Ej.  Now  that  the  are  known  as  function  of  Ej,  it  is 
relatively  easy  to  show,  using  Equation  (B.l),  that  the  second  deriva- 
tive of  W evaluated  at  equilibrium  (e3  = 1)  is 


32w  3 £ 2 

~ =i7-  (ctl-c12) 

0 


Having  defined  a deformation  Ej  that  is  related  to  the  elastic 
constants  Cjj  and  C12,  we  now  consider  a deformation  £2  which  can  be 
related  to  the  C..  elastic  constant.  The  constant  C, , is  related  to 

*♦4  4 4 

a shear  that  changes  the  angle  between  the  z and  x axis.  An  alterna- 
tive way  of  viewing  this  shear  is  that  vectors  along  (1,1,1)  expand 
while  vectors  that  are  perpendicular  contract  or  vice-versa  (see 
Figure  18)*  Hence,  the  deformation  e , takes  the  vectors 


I 
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Sj  = (1,1,-D,  s2  = (1,-1,0),  s3  = (1 ,0,-1) 
and  transforms  them  to 

sj  = A (1 , 1 ,-l) , sj  = B (1 , -1 ,0) , sj  = B ( 1 , 0 , -1 ) 


(B.8) 


(B . 9) 


Taking  Equations  (B.3),  (B.4) , (B.8),  and  (B.9)  't  may  be  shown  that 


A + 2B  - 3 


and 


e„  = e5  = es  = T (A_B) 


(B.10) 


_i /2 

It  may  further  be  shown  that  B = A in  order  to  ensure  volume 

conservation  of  the  primitive  unit  cell.  Finally,  if  A is  replaced  by 
i 

e2  , then 


del 


3 l2C, 


(B.ll) 


where  subscript  0 indicates  the  derivative  is  evaluated  at  equilibrium. 


We  have  now  defined  deformations  £j  and  e2  that  are  related  to 
Cjj  - C12  and  C44  through  the  second  derivative  of  the  elastic  energy 
density.  To  evaluate  these  constants, it  is  necessary  to  derive  an 
expression  for  W,  which  is  the  energy  per  ion.  This  energy  is  com- 
posed of  two  parts  that  will  give  nonzero  contributions  to  the  deriva- 
tives  in  EXjuations  (B.7)  and  (B.ll).  W consists  of  W , the  conduction 

6 S 0 S 

electron  energy,  and  W , the  electrostatic  energy  per  ion.  W 


represents  the  average  electrostatic  energy  per  ion  in  a crystal, 
where  the  crystal  is  visualized  as  a regular  array  of  point  positive 
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es 

charges  immersed  in  a compensating  uniform  electron  sea.  W will  be 

evaluated  with  the  aid  of  a technique  developed  by  Ewald  for  producing 

ce 

convergant  lattice  sums.  The  other  contribution,  W , is  the  struc- 
turally dependent  part  of  the  conduction-electron  energy.  Actually, 

Cg 

W has  been  encountered  previously  as  the  quantity  E in  Equations 
(18)  and  (A. 19).  It  is  the  second  order  contribution  to  the  conduction 
electron  energy.  The  zero  and  first  order  contributions  to  the  elec- 
tronic energy,  being  structurally  independent,  do  not  contribute  to  the 

GS  CG 

derivatives  in  Equations  (B.7)  and  (B.ll) . The  sum  of  W and  W give 
the  energy  of  a primitive  unit  cell  in  the  crystal.  The  square  of  the 
structure  factor,  S(q),  in  Equation  (18)  is  responsible  for  the  ex- 

CG 

plicit  dependence  of  W upon  the  crystal  structure.  For  a perfect 
lattice,  the  absolute  square  of  S (q)  is 


| S (q)  | 2 = — £ e-iq*X* 

N2  l 
= oq,G 


(B.12) 


In  Equation  (B.13),  x^  represents  a lattice  vector,  G a reciprocal 
lattice  vector  and  N is  the  number  of  unit  cells  in  the  crystal.  Hence 
Equation  (18)  becomes,  with  the  aid  of  Equation  (B.12) 


WCe  = l F(G)  . 
-+ 

G 


(B.13) 


es 

The  electrostatic  energy  W is  a rather  complicated  expression 
so  we  shall  give  some  of  the  details  of  its  derivation.  The  charge 


density  of  the  crystal  can  be  written  as 


* 

z e 


l e"iq*X*  - 
a 


6q,0 


(B. 16) 


The  Coulomb  potential  associated  with  p(x)  is 


c 4TT  iq*x 

|(x)  =2.  — P(q)  e 
q q2 


* , iG*x 

4tt  z e y e 

Z „2 


(B. 16) 


where  is  the  primitive  cell  volume,  z*e  the  charge  of  an  ion  in  the 
crystal,  and  the  prime  on  the  G sum  denotes  that  the  G = 0 term  is 
omitted.  However,  one  should  note  that  Equation  (B.16)  does  not  repre- 
sent a converging  sum.  We  must  resort  to  a method  due  to  Ewald  to 
obtain  a convergent  sum  equivalent  to  Equation  (B.16) . One  may  easily 
verify  that 


♦* 


ac 


(B.17) 


Equation  (B.17)  is  identical  to  Equation  (B.16) . Factoring  the 
integral  into  different  domains  of  integration 


x 4tt  z e j 

*(x}  m ~nT~)\ 


' e'^  + iG‘X  d£  ♦ 


.“° 5 + iG*x  , 


(B.18) 


Equation  (B.19)  represents  the  Coulomb  potent  i 1 at  x due  to  all  chat 


in  t ho  crystal.  Subti acting  the  potential  due  to  the  point  charge  at 


♦ 

from  Equation  (It.  19),  letting  tl 
obt a 1 n 


1 . » * 

, and  sett  1 ng  x - x,,  wo 

2/n 


4'’Ue, ) 


, z*e  et fc  n 
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i*i 


•1  nz  * o 

ihr’ 


which  is  tlie  jvtential  at  x ^ - tl  due  to  all  charges  in  the  crystal 
except  the  point  ion  at  x’  « 0. 

In  the  abovq  wo  determined  the  }\itantial  4'  which  was  the  poten- 
tial  due  to  all  but  one  charge  in  the  crystal . The  single  charge  is 
omitted  in  deriving  t'quot  ion  111. 201  since  it  is  not  correct  to  have  a 
charge  interact  with  itself  when  the  interaction  energy  ot  the  system 
is  computed.  ^ may  be  lelated  to  W as  follows 


C 6 S 

Expressions  ( B . 13)  for  W , and  (B.21)  for  W can  be  put  into 


either  liquation  (B.7)  or  (B.ll)  to  compute  the  elastic  constants 


9*W  3V'e  3JWeS 
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-G"  -2-  + -L- 
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lB.22) 


G' , G",  x£,  x£'  are, respectively,  the  first  and  second  derivatives  of 


the  reciprocal  and  real  lattice  vectors  with  respect  to  Cj  or  C?  and 


the  other  quantities  have  been  previously  defined.  With  Equation 


(B.22),  we  have  all  the  information  needed  to  calculate  the  elastic 


constants  C,,  - Cw  and  CU1)  for  a bcc  crystal.  The  quantities  G', 


G",  x£,  x"  can  be  evaluated  easily  since  x£(t\)  and  G(e  J are  known  be- 


cause the  strain  matrix  e has  previously  been  evaluated  for  the 


deformations  t , and  c,. . 

I 4 
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APPENDIX  C 

THE  ORTHOGONALIZATION  HOLE 

The  orthogonal izat ion  hole  (OH)  density  is  one  component  of  the 

conduction  electron  charge  density.  In  this  appendix  an  exact  expres- 
OH 

sion  for  n (r)  is  derived  and  compared  with  the  standard  approximations 

OH  “*■ 

of  Harrison  (1966).  The  core  shift  due  to  the  interaction  of  n (r) 

with  the  core-electron  density,  V°^w  is  computed  below  with  both  the 

OH 

approximate  and  exact  forms  of  n (r)  . 

The  conduction-electron  charge  density,  nconc3(r)  , is  obtained  by 

c 

talking  the  absolute  square  of  ip  (r)  and  summing  over  all  occupied  k 

k 

states. 


cond 
n (r) 


l 4> * ( r ) ( 1 - P*)  (1  - P)^(r) 

+ k k 

k 


(C.l) 


Equation  (7)  was  used  to  obtain  expression  (C.l)  for  the  conduction- 

electron  density.  Substituting  equation  (6)  into  (C.l)  and  keeping 

cond 


only  zero  or  first  order  terms/  n 


cond 


£>  = l lb  - 


-ik»  r 


P k> 


(r)  becomes, 


ik*r 

e --‘•i  * 


AT 


<k I P*  + <k|P*,p|k> 


1 c ( * 

+ n i 

q*> 


.r*  +iq*r  ,7*.  -iq*r/ 

a (k)  e + a (k)  e 


(C.  2) 


Plane  wave  states  not  being  operated  upon  by  the  projection  operator 
have  been  shown  explicitly. 


Wr 


The  fact  that  a (k)  » 1 has  also  been  used  in  the  derivation  of  Equation 
q-0 

OH  ■* 

(C.2),  the  2nd,  3rd,  and  4th  terms  of  Equation  (C.2)  comprise  n (r) 

sc  -*■ 

while  the  last  term  constitutes  the  screening  charge  density  n (r) . 

The  OH  density  is  essentially  determined  once  the  core-electron  wave 

SC 

functions  have  been  specified,  whereas  n (r)  represents  the  self- 
consistent  response  of  the  electron  gas  to  the  perturbing  potential  W(r) . 

It  is  convenient  to  make  the  following  substitutions  in  Equation 

(C.2).  The  core  eigenstates  in  the  projection  operator  P = \ |nX><n£| 

nX 

are  explicitly 


|nX>  = P^U)  Yt(0,«l»/r  . (C.3) 

Also,  the  plane  waves  may  be  expanded  as 

_ -*  -► 

elk*r  - l 1 4 n (2X+1) ] '/j  j . (kr)  ilY®(0,4>)  , (C.4) 

X t x. 

where  P „ (r)  is  the  radial  part  of  the  core  wave  function,  j.  is  a 
nX  x 

spherical  Bessel  function,  and  Y^1  is  a spherical  harmonic.  With 

OH  -*■ 

Equations  (C. 3)  and  (C.4)  in  Equation  (C.2),  n (r)  can  be  expressed  as 


OH  - "4n» 

n (r)  = 


(24+l)<k  ni>P  „(r) 
nX 


( 2tt ) 3 I o nX 


,17  , Pn't<rl  , 

^ h <kr>  - l.  k dk 

n 


The  quantity  <k | n£>  is  known  as  an  orthogonality  coefficient  and  is 
given  by 
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<k  n£> 


4n  f . 

^0  J ^ 


j^(kr)  P^tr)  rdr 


cond  OH 

The  expressions  (C.2)  for  n (r)  and  (C.5)  for  n (r)  are  not 

correct  as  they  stand  because  the  conduction  wave  functions  are  not 
cond 

normalized.  If  n (r)  is  integrated  over  all  space,  one  finds  for 
the  total  charge  of  a primitive  unit  cell 


nCOndU)  d3r  = 1 - <k|p|k> 


Hence,  n^H(r)  and  n'’°nd(r)  should  be  divided  by  the  factor  1 - <k|p|k> 
to  ensure  the  correct  normalization  for  the  conduction  charge  density. 
When  the  normalization  factor  in  Equation  (C.7)  is  included  in  Equations 
(C.2)  and  (C.5) , one  finds  that  there  are  effectively  Z*  electrons  in 
nSC  (r)  and  Z - Z*  electrons  in  n0H(r),  where 


<k  P k> 


Z is  the  effective  valence  of  the  metal,  and  Z^  is  the  classical 
valence  (e.g.,  = 1 for  Li,  Na,  K) . The  effective  valence  Z*  will 

always  be  slightly  greater  than  Z since  <k|P|k>  is  a positive  definite 
quantity. 

OH  -*■ 

The  charge  density  n (r)  (properly  normalized)  will  give  rise  to 

OH 

two  effects.  Firstly,  its  Coulomb  potential  V (r)  will  contribute  to  the 

c -*■  OH  -*■ 

crystalline  potential  V (r) , and  secondly,  n (r)  will  interact  with  the 
core 

core  charge  n (r)  to  shift  the  core  eigenvalues  by  an  amount  called 
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opw  OH 

V . The  potential  V (r)  was  evaluated  by  performing  a standard 

OH  •* 

Coulomb  integral  over  the  spherically  symetric  charae  density  n (r) . 
Then  the  energy  shift  of  the  n , 9. ' th  core  state  is  oiven  by 


v°pw  = e2 
n«. 


(r)  V°H(r)  dr 


The  matrix  element  of  V (r)  is  also  needed  in  the  evaluation  of  the 
pseudopotential,  and  with  the  aid  of  Poisson's  equation  can  be  shown 


<k  + q|v°H(r) |k>  = ~ — n°H(q)  , 

q2fl 

1 0 


(C.10) 


where  n0,‘ (q)  is  the  Fourier  transform  of  n°H(r) 


Equation  (C.5)  constitutes  an  exact  expression  for  n (r)  when  the 
normalization  factor  (C.7)  is  included.  All  calculations  in  this 


thesis  were  done  with  the  exact  form  of  n (r) . Harrison  (1966) 
originally  proposed  using  approximate  treatments  of  n°H(r)  and  V°pw 


because  of  the  difficulty  in  evaluating  the  exact  n (r)  . However  as 
will  be  shown  these  approximations  are  rather  severe  and  can  lead  to 


significant  errors.  Harrison  used  the  following  approximation  for 

OH 

the  matrix  element  of  V (r) 


core 

+ -*■  | OH , r*  H 4ne2  (q  * 

<k+q  V (r)  |k>  (Z  - 7,  ) . 

' core  V 

q “0  n (q=0) 


(C.ll) 


Equation  (C.ll)  is  based  on  the  assumption  that  the  OH  is  distributed 
in  the  same  manner  as  the  core-electron  charge  density.  The  OH 


9 


core  shift  on  the  other  hand  was  approximated  by 


yOpw-H 


*o  v 

} 2(24+1) 

tr  irC 


fP' (r)dr 

n_e 

t 


(/  <k|n4>2  k*'dkl 


(C . 12) 


In  Equation  (C.12),  the  OH  is  approximated  as  a point  charge  of  magnitude 

Z “ located  at  the  nucleus.  Clearly,  this  is  a rather  crude 
V 

approx i mat  ion . 

It  is  instructive  to  compare  the  exact  treatment  of  the  OH  [i.e., 
Equations  (A. 5),  (C.9),  (C .10))  with  the  approximate  treatments  li.e., 

Equations  (C.ll)  and  (C.12)).  Table  i>  shows  calculated  exact  and 

opw  opw 

approximate  values  of  V for  Li,  Na , and  K.  The  approximate  V 

is  roughly  a factor  of  2 to  3 too  large  and  is  not  a function  of  n and 

4 as  it  should  be.  Cutler  et  al.  (1975)  showed  for  zinc  that  the 

Ol'W 

differences  in  these  two  treatments  of  V can  lead  to  significant 

changes  in  calculated  lattice  dynamical  properties  such  as  the  phonon 

opw 

spectra.  Li,  Na , and  K also  are  affected  by  t tie  differences  in  V 

OH  » 

The  approx  invite  and  exact  treatments  of  n (r)  yield  quite  different 

OH  * 

results.  Figure  19  shows  the  exact  n (r)  for  Li,  Na , and  K that  were 
OH  * OH 

used  to  calculate  V (t).  Notice  that  n (r)  changes  sign  in  certain 

OH  OH  * 

regions  of  r.  When  the  approximate  V (r)  is  calculated,  n (r)  is 

assumed  to  have  the  same  spatial  distribution  as  the  core  electrons. 

OH  * 

Hence,  the  approximate  n (r)  will  not  exhibit  any  changes  in  sign  and 
is  therefore  not  even  in  qualitative  agt eement  with  the  exact  treatment. 


k 


Table  6. 

Exact  versus 

approximate  V°pw  for 

Li,  Na,  and  K. 

Exact 

Approximate 

Energies  in 

v°pw 

V°PW 

Rydbergs 

Li 

-0.0743ls 

-0.2659 

Na 

-0. 0594 ls 
-0.08682s 
-0.08512P 

-0.1870 

K 

-0.0988ls 
-0. 11172s 
-0.11032P 
-0 . 12353s 
-0.12323p 

-0.2029 
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APPENDIX  D 

SCREENING  OK  THE  PSEUDOPOTENTIAL 


This  appendix  presents  a brief  derivation  of  the  screening 
contribution  to  the  pseudopotential.  The  conduction-electron  charge 
density  [i.e.,  Appendix  C,  Equation  (C.2)l  consists  of  the 
orthogonal ization  hole  density  which  was  discussed  in  detail  in 
Appendix  C and  the  screening  charge  density  which  is  explicitly 


sc 

n 


<r> 


+ l fa*  HO 
q*0l  q 


• r 


-t 


a (it) 

q 


(D.  1) 


► cond 

The  coefficients  a (k)  were  previously  used  to  expand  \p  (**'  and 

a k 

A (r)  in  terms  of  opw 1 s and  plane  waves,  respect ively . First  order 

k 

perturbation  theory  applied  to  Equation  (4),  the  pseudopotential 

-> 

equation,  yields  a^(k)  as 

q 


a^  (k) 

q 


<k+q  |W|  k> 

— (k‘ - | k+q | * ) 


Equation  (D. 1)  can  be  rewritten 


(D.  2) 


sc 

n 


(r) 


(2n)  3 


d ‘k 


l 

q^O 


a (k) 
_-q 


a (k) 

q _ 


(D.  3) 


where  the  sum  over  k has  been  converted  to  an 
standard  manner.  Now,  if  one  notes  that  av(k) 

q 


inteoration  in  the 

* 4 

= a + (-k)  and  that 

-q 


* * * "*■ 

a ^(k)  may  be  replaced  by  a due  to  the  symmetric  k integration, 


n^c(r)  may  be  conveniently  expressed  as 


sc,  \ 
n (r)  = 


(27!)  3 


i + 2y  


q^O  Jf\ 


^(k2-|k+q|  2) 


Interchanging  the  orders  of  the  integration  and  summation  and  express- 
sc  * 

ing  n (r)  as  a Fourier  transform  leads  to 


sc  . 
n ( r) 


= l nSC(q) 


sc,-*,  in*r 

(q)  e ' 


whe  re 


r 3 > -*  | , 

d k <k+q  W kN 


(2-'  3 A2  2 , 


q / 0 


— (k2- ! k+q I 2 ) 


SC  ► 

n (q)  = 


>tt)  3 


[■  d3k  = Z q = 0 


When  the  correct  normalization  of  the  conduction  wavefunction  is  used 
(see  Appendix  C),the  q = 0 Fourier  coefficient  will  have  the  value 
Z /i.  . The  Fourier  components  of  n (r)  for  nonzero  q represent  the 

perturbative  response  of  the  conduction-electron  qas  to  the  perturbing 
potential  W(r)  . We  shall  now  describe  how  the  response  of  the 
conduction  gas  to  the  perturbing  potential  W(r)  can  be  made  self- 
consistent.  We  will  also  show  how  exchange  and  correlation  effects 
between  conduction  electrons  can  be  incorporated  into  the  description 


of  screening. 


I 
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In  Equation  (D.6),  it  is  possible  to  make  the  screeninq  self- 
consistent.  by  usinq  Poissons  equation  and  by  separating  W(r)  into 

a bare  unscreened  potential  v/'(r)  and  the  Coulomb  (Potential  of  the 

SC  sc 

screeninq  electrons  w'  (r) . The  Fourier  transforms  of  n (r)  and 

sc 

W (r)  are  related  by  Poissons  equation. 


w’1  (q)  = nbC(q)[l  - G(q)] 


(D.7) 


The  Coulomb  interaction  lias  been  modified  by  the  exchange  and  correla- 
tion interaction  between  the  electrons  in  the  conduction  qas  and  this 
is  explicitly  taken  into  account  by  the  extra  factor  of  1 - G(q)  in 

Equation  (D.7)  (Singwi  et  a.1 . , 1968) . Using  Equation  (D.7)  in  Equa- 

b sc 

turn  (D.6)  and  separatino  w into  W + w , one  obtains 


VJSC(q)q2  tl-G(q)  1 


4H  e 


4 sc.  . 

= W (q) 

(211)  3 


d3k 


2m 


(kMkVqV) 


(2n)  3 


<lt-tq|wb|k>  d3k 


( D.  8) 


In  Equation  (D.8),  w (<i)  is  factored  out  of  the  integrand  because 
-► 

it  .s  independent  of  k.  Takinq  note  of  the  fact  that  (Bardeei\  1937) 


r 


d3k 


f-(K?- 


(1  - e (q)  ] 


(D.9) 


I k+q | 2 ) 


who  i o 


,2  ^ i -q  /4  2k  *1] 

t (,])  - 1 ♦ k t 1 I’n 

-xfV  ‘ 11  -S  'l 


(0. 10) 


one  may  show  th.it 


80  , .«e;  lWU''H  d’k 

W (q)  . 

.'m 


(n.  in 


The  tunot  ion  i (q)  in  t ho  nt.it  io  di  <•  loot  t i o rosixinso  function  toi  .in 
int e i .not  1 nq  election  qas.  It  in  related  to  the  ttee  elect  toil  v’l  no- 
o.i  I l«nt  It. 1 1 1 t ee  n t ,i  t i o .1  i e I eot  i i o t unot i on  i (q ) hy 


(q)  - 1 l l (»(q ) 1 t I (q)  ~ 1 I 


(H.  1 .') 


Kqu.it  ion  (0.  11)  toi  tho  sevooninq  oont  v iluit  i on  to  the  psoudo- 

ixitenti.il  indicates  tli.it  one  may  express  the  Hcieeiuiui  potenti.il  in 

terms  ol  the  Kuo  peiturliinq  (xitenti.il  ami  a dioloctiic  nerooniiui 

function  i *(q).  The  matrix  ('lenient  i't  the  h.ne  oi  mmol  eeni'il  |'.it  t ot 

► » | h.' 

the  (':ieiulo('otent  i al  'k*.(|W  |k-  in  oivon  in  pqu.ition  ill'.  In  local 

h 

pn«'ll(lopot  ('lit  l al  cali'iil  at  ions,  t he  in.it  i l x element  ol  W in  tactoiod  out 

i't  the  i lit  oq  i and  in  Pqu.ition  (IV  11)  and  a veiy  simple  iet.it  i»»n  i •» 

h no 

thus  obtained  between  W lq)  and  W (q) . 


.... 

W (q)  ~ - w (.() 


(!'.  1 1) 
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Although  this  is  a widely  used  approximation,  it  was  not  used 
in  any  calculations  in  this  thesis  because  the  bare  potential  matrix 
elements  are  not  independent  of  k as  implied  by  this  approximation. 

C -> 

The  matrix  element  of  V (r)  in  Equation  (14)  is  independent  of  k and 
could  be  factored  out  of  the  integrand  in  Equation  (D.ll),  but  the 
remaining  terms  in  Equation  (14)  are  nonlocal ly  dependent  on  both  the 
magnitude  and  direction  of  k.  The  local  approximation  siqni f icantly 
alters  the  screening  field  and  consequently  affects  the  predictions  of 
lattice  dynamical  properties.  Therefore,  the  full  nonlocal  form  of  the 
pseudopotential  was  used  in  all  our  calculations. 
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APPENDIX  E 

DERIVATION  OF  THE  LIQUID  METAL  RESISTIVITY  EXPRESSION 

Ziraan's  formula  (Ziman,  1961)  for  the  electrical  resistivity  0 of 

liquid  metal  was  given  in  Eauation  (39)  where  k , v , n , w(q)  , 

f f o 

S (q)  , e,  and  q are  the  Fermi  wavenumber,  Fermi  velocity,  avaraae 

number  density,  pseudopotential  form  factor,  liquid  structure  factor, 
electronic  charae,  and  momentum  transferred  in  a scattering  process/ 
respectively.  We  will  now  briefly  discuss  Equation  (39)  and  indicate 
what  approximations  are  inherent  in  this  formula. 

The  current  density  J for  a system  of  charged  particles  in  a 
volume  subject  to  an  electric  field  F along  the  z axis  can  be 
expressed  in  terms  of  the  distribution  function  f(k). 


= - y 

■ o L 


e V f (k) 


(E.l) 


Converting  the  sum  over  occupied  k states  to  an  integral  in  the 
standard  manner  yields, 

V f (k)  d3k  . (E.2) 

z 


If  the  distribution  function  is  a linear  function  of  F,  then  the  scalar 
conductivity  0 and,  subsequently,  the  resistivity  p can  be  determined 
from  the  relation  J = OF.  To  find  J or  0,  we  must  first  calculate 
f (k)  . For  a system  that  is  in  a steady  state,  Liouville's  equation 
states  that. 


dfOO/dt  = 0 


(E.  3) 
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It  is  convenient  to  reexpress  the  total  time  derivative  as  a sum  of 
partial  time  derivatives.  Ignoring  diffusion  and  assuming  zero 
magnetic  field,  the  distribution  function  f(k)  will  change  with  time 
because  of  the  field  F and  scattering  processes.  Hence,  Equation 
(E.3)  is  equivalent  to 


3f (k)  3f  (k) 
3t_.  ,4  3t  - 

field  scat 


(E.4) 


where 


3f  (k) 

3t  . 

f leld 


3f (k)  kz  eF 
3k  k h 


and 


3f  (k) 


3t 


scat 


/f(k')  [l-f(k)]P  ,,d V 


/f (k) [1-f  <k')  JP  ,d  V 


(E.5) 


(E.6) 


The  function  Pkk>  is  the  probability  for  a transition  from  the  state  k 
to  the  state  k'.  Combining  Equations  (E.4),  (E.5),  (E.6)  and  writing 
f(k)  as  the  sum  of  an  equilibrium  contribution,  f°  (k) , and  a contribution 
g(k)  due  to  the  electric  field  results  in 

f - /hid  - ,(k)iPkk,d'k'  . (e.7, 

"V 

Equation  (E.7)  was  obtained  by  keeping  only  terms  linear  in  F (note 
that  g(k)  is  proportional  to  F so  that  we  are  limited  to  small 
signals]  and  assuming  that  pkk  > = 


The  electric  field  changes  the  component  of  an  electron's  k 

vector  by  an  amount  eFt/h  in  time  t.  However,  the  scattering  processes 

tend  to  offset  the  momentum  gain  of  the  electron  caused  by  the  electric 

*♦ 

field.  If  the  electric  field  were  turned  off,  we  would  expect  J to 
decay  exponentially  with  a characteristic  relaxation  time  T that  is 
determined  by  the  scattering  mechanism.  Similarly,  the  Fermi  distri- 
bution would  change  in  time  due  to  the  electric  field  if  the  scattering 
processes  did  not  stop  the  electrons  momentum  from  increasing  indefinitely. 
It  is  easily  seen  that  the  perturbation  part  of  the  distribution  function, 
g(k),  can  be  expressed  as 


3f°(k)  eF  z _ 

g(k)  3 -W~TTX 


where  T is  a characteristic  relaxation  time  of  the  system.  To  obtain 
a solution  for  g(k),  we  assume  the  Born  approximation  for  P , 

KK 


7 Kk'l!SS(W>  ' 


where  W^»  is  the  matrix  element  of  the  scattering  potential  between 
the  states  k and  k'.  Using  Equations  (E.8)  and  (E.9)  in  Equation  (E.7) 
yields  the  following  result  for  the  relaxation  time  T 


I*  milk ' 


d cos  P 


(E.10) 


where  0 is  the  angle  between  k and  k'.  The  use  of  the  Born  approximation 
means  that  k and  k both  lie  on  the  same  energy  sphere  (i.e.,  |k|  *=  |k  |). 
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The  expression  for  I qives  us  complete  knowledqe  of  f (k ) and  it 
is  now  a trivial  matter  to  formally  calculate  the  current  density  J 


, ( v k 
2e  2 2 


2 2 eP  H f ( k ) i 

T J kf  5?  T*  1 d k 


(E.ll) 


(2H) 

Assuminq  that  3f°/3k  is  a sharply  peaked  function  of  the  Fermi  surface 


4e2  k’  T 

J = — F 

" m ( 2tt)  2 3 


(E.12) 


It  is  obvious  from  Equation  (E.12)  that  the  resistivity  p for  the 
metal  is 


P = 


3(2n)‘  m 
4e2  kj  T 


(E. 13) 


which,  after  some  alqebra,  coordinate  chanqes,  and  factori2inq  of  the 
pseudopotential  matrix  element,  in  the  expression  for  1/T  reduces  to 
Equation  (39). 

Some  qeneral  comments  may  be  in  order  at  this  point.  In  the 
derivation  for  p, it  was  assumed  that  there  was  a sharp  Fermi  surface 
so  that  3f°/3t  was  a delta  function  at  this  surface.  However* one  can 
estimate  that,  if  the  temperature  blurrinq  of  fp  was  taken  into  account, 
the  correction  to  p would  be  rouqhly  four  orders  of  maqnitude  smaller 
than  the  results  obtained  with  the  sharp  Fermi  surface  assumption. 
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